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SUMMARY 


In  the  present  context,  a  compliant  robotic  structure  is  one  or  more 
continually  flexible  beams  whose  motion  can  be  controlled  to  manipulate 
objects.  Typically,  a  single  beam  is  comprised  of  many  hollow  tube  elements 
placed  end-to-end,  where  each  tube  element  moves  in  a  particular  way  when 
pressurized  with  air.  Each  tube  element,  made  of  a  reinforced  polymeric 
material,  is  designed  with  directional  stiffness  properties  so  that  tube 
motion  (bending,  twisting,  extension,  coiling)  is  predictable  for  a  given 
external  loading  and  internal  pressure  history.  The  basic  objectives  are  to 
develop  theories  for  the  motion  and  control  of  selected  single  elements  and 
multiple  element  systems  (beams)  to  be  used  a  light-weight,  fast-acting 
manipulators. 

To  meet  the  objectives,  the  first  year  of  this  three-year  project  was 
focused  on  three  separate  yet  related  topics.  The  first  involved  studies  of 
an  elephant's  trunk  from  the  viewpoint  of  a  compliant,  continuously  flexible 
manipulator.  The  relationships  among  the  load-carrying  capacity,  material  (or 
muscle)  properties  and  trunk  geometry  were  investigated  and  are  reported  in 
Chapter  I.  Using  clues  from  this  animal  study,  particularly  the  observations 
of  how  a  trunk  coils  around  an  object  and  then  lifts  it,  the  second  study 
involving  the  nonlinear  mechanics  of  bending  and  coiling  tube-type  element 
systems  was  carried  out.  Those  results  are  reported  in  Chapter  II.  The  third 
study  employed  the  linear  theory  of  elasticity  to  predict  the  load  and 

deformation  behavior  of  cylindrical  tubes  with  orthotrophy  in  the  form  of 

c 

constant  angle  helves.  The  quantatlve  results  given  in  Chapter  III  show  that 
such  a  tube  element  can  be  used  as  a  torsion  actuator,  or  a  manipulator  that 
twists  about  its  longitudinal  axis  when  pressurized. 
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I.  MECHANICS  OF  AN  ELEPHANT  TRUNK: 
THEORY  AND  MEASUREMENTS 


ABSTRACT 

Elephant  trunks  are  versatile  and  strong.  Their  mechanical  design  may 
therefore  be  used  as  a  source  of  ideas  for  the  design  of  versatile,  strong 
robotic  manipulators.  In  this  paper  a  mathematical  model  of  an  elephant  trunk 
lifting  a  weight  is  developed  based  on  assumptions  parallel  to  those  used  in 
the  analysis  of  prestressed,  linear,  composite  beams.  Data  on  overall  trunk 
geometry  obtained  during  weight-1 ifting  and  data  for  the  trunk  muscle  distri¬ 
bution  at  a  cross  section  are  used  to  calculate  an  apparent  tangent  modulus 
for  the  trunk  tissue.  During  a  change  in  trunk  curvature  for  which  the  uni¬ 
form  component  of  longitudinal  prestrain  remains  a  constant  30  percent,  a 
typical  value  for  the  apparent  modulus  is  of  the  order  of  lO^N/m^.  This  study 
represents  an  initial  step  toward  more  refined  analyses  that  rationally  relate 
the  loading,  the  overall  trunk  geometry  and  the  muscle  structure  while  satis¬ 
fying  equilibrium  conditions. 


INTRODUCTION 


The  trunk  of  an  elephant  is  a  tapering  muscular  structure,  completely 
lacking  rigid  materials,  of  nearly  circular  cross-section.  It  is  versatile 
and  strong;  trunks  can  shorten,  extend,  bend,  and  twist,  and  elephants  have 
been  trained  to  lift  and  manipulate  masses  of  up  to  300  kg  with  their  trunks. 
This  study  is  part  of  a  larger  study  exploring  the  mechanics  of  trunk  motion 
and  strength,  and  aimed  at  discovering  principles  of  trunk  design  that  may  be 
used  in  the  design  of  similar,  continuously  flexible,  versatile,  strong 
robotic  manipulators.  A  further  discussion  of  animal  muscular  hydrostats  and 
robotic  mechanics  is  given  by  Wilson  (1984). 

The  purpose  of  the  present  study  is  to  correlate  original  data  for  the 
overall  geometry  of  a  trunk  lifting  a  load  with  the  reactive  forces  of  the 
trunk  muscles,  to  measure  the  average  longitudinal  strain  of  a  trunk  lifting  a 
weight,  and  then  to  calculate  an  apparent  modulus  of  trunk  tissue.  The 
apparent  modulus  is  defined  as  the  average  tangent  modulus  for  all  of  the 
longitudinal  and  oblique  muscle  tissue  at  a  given  cross  section.  This  average 
tangent  modulus  Ea  represents  the  slope  of  the  stress-strain  curve  in  a  uni¬ 
axial  test  in  which  this  whole  muscle  group  has  a  uniform  longitudinal  pre¬ 
strain  cv.  The  method  used  to  calculate  the  modulus  is  one  that  does  not 
require  dissection  of  the  animal  involved,  and  may  be  used  by  zoologists 
Interested  in  the  apparent  modulus  of  other  animal  appendages.  This  modulus 
calculated  for  the  trunk  may  be  used  as  a  guide  in  choosing  a  material  for 
trunk-like  robotic  manipulators. 

The  size  and  distribution  of  the  trunk's  four  basic  muscle  masses— the 
radial,  the  longitudinal,  and  the  two  oblique  layers— probably  account  for  the 
elephant's  ability  to  grasp  and  manipulate  loads  with  Its  trunk.  Although  the 
role  of  each  muscle  type  in  trunk  movement  and  load  manipulation  is  not  fully 


understood,  it  is  probable  that  one  can  draw  parallels  to  the  mechanical 
behavior  of  squid  tentacles  and  vertebrate  tongues  as  outlined  by  Kier  (1982). 
In  this  case  contraction  of  selected  oblique  muscles  would  lead  to  torsion  of 
the  trunk  about  the  longitudinal  axis,  uniform  contraction  of  radial  muscles 
would  lead  to  uniform  trunk  lengthening,  and  contraction  of  longitudinal  and 
radial  muscles,  acting  separately  or  together  on  one  side  of  the  trunk,  would 
lead  to  trunk  bending. 

In  this  analysis  the  longitudinal  trunk  loads  (self-weight  and  payload) 
are  assumed  to  be  supported  solely  by  the  tensions  of  the  longitudinal  and 
oblique  muscles.  The  corresponding  longitudinal  strain  has  two  components:  a 
relatively  larger,  uniform,  contraction  prestrain  (observed  to  be  20  to  35 
percent  before  and  during  lifting);  and  a  relatively  small  strain  (two  percent 
maximum)  due  to  trunk  bending  and  curvature— as  the  payload  is  lifted.  In 
achieving  the  uniform  contraction  prestrain,  the  stress-strain  behavior  is 
probably  very  nonlinear.  The  tangent  modulus  (the  slope  of  the  stress-strain 
curve)  may  Increase  dramatically  since  muscles  are  much  stiffer  after  con¬ 
traction  (Yamada  and  Evans,  1970).  For  the  small  strains  associated  only  with 
trunk  bending  and  small  changes  in  trunk  curvature  between  two  states  of  lift¬ 
ing  the  payload,  the  analysis  is  based  on  classical  beam  theory.  That  is,  the 
trunk  is  assumed  to  be  composed  of  a  linear,  composite  material  (the  longi¬ 
tudinal  and  oblique  muscles)  for  which  the  apparent  modulus  on  the  side  of  the 
neutral  axis  under  the  least  muscle  tension  is  the  same  as  the  apparent 
modulus  on  the  other  side  under  the  most  muscle  tension.  The  apparent  modulus 
of  this  analysis,  then,  represents  the  average  value  of  the  tangent  modulus 
for  the  muscles  Involved  In  Incremental  bending,  given  a  uniform  prestrain. 

The  analysis  implies  that  during  Incremental  bending  there  is  no  change  in 
chemical  structure  of  the  material  with  load  other  than  the  strains  imposed  on 


the  chemical  bonds  by  the  load.  Clearly,  the  trunk  of  a  live  elephant  may  not 
exactly  satisfy  all  of  the  assumptions  and  conditions.  However,  local  varia¬ 
tions  in  the  modulus  of  the  components  of  muscle  cells  and  fibrous  connective 
tissues,  and  other  local  variations  such  as  possible  shifts  in  the  neutral 
bending  axis  and  possible  nonlinear  bending  strain  distributions,  would  have 
little  effect  on  either  the  average  apparent  modulus  for  the  load-carrying 
tissue  acting  as  a  group,  or  on  the  overall  mechanical  performance  of  the 
trunk.  The  approach  herein  represents  the  first  step  towards  a  more  refined 
and  accurate  analysis  of  the  trunk. 

Briefly,  the  analysis  and  measurements  proceed  as  follows.  As  described 
in  the  Mathematical  Model  and  Analysis  section,  a  linear  strain  distribution 
at  a  cross  section  and  force  equilibrium  were  used  to  calculate  the  location 
of  the  neutral  axis  of  the  trunk.  To  do  this  the  cross-sectional  areas  of 
muscles  were  obtained  from  a  drawing  of  a  trunk  cross-section  and  a  discussion 
of  detailed  trunk  morphology  presented  in  the  classical  study  of  the 
elephant's  head  by  Boas  and  Paulli  (1908).  The  drawing  of  the  cross-section 
was  then  estimated  to  coincide  with  a  particular  location  along  the  length  of 
the  trunk.  At  this  location  the  radii  of  curvature  po  and  P}  at  two  lift 
positions  were  measured  from  films  we  took  of  an  elephant  lifting  a  known 
weight  with  its  trunk.  The  second  moment  of  area  of  the  muscles,  I,  was 
determined  from  the  size  and  distribution  of  the  muscles  in  the  Boas  and 
Paulli  cross-section.  Based  on  moment  equilibrium,  pq»  pj,  and  I  were  then 
used,  as  described  below,  to  calculate  the  apparent  modulus  of  trunk  tissue. 

This  paper  is  divided  into  four  parts.  In  the  first  two,  the  Muscle 
Morphology  section  and  the  Lifting  Experiments  section,  we  describe  the 
sources  of  our  measured  data.  We  develop  our  mechanical  model  of  the  trunk  in 
the  Mathematical  Model  and  Analysis  section.  We  then  present  our  calculated 
data  in  the  Numerical  Results  and  Discussion  section. 
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MUSCLE  MORPHOLOGY 


Figure  1  presents  a  schematic  drawing  adapted  from  Boas  and  Paul  1 i  (1908) 
showing  the  placement  and  orientation  of  the  four  basic  muscle  masses  and  the 
two  main  tendinous  masses  in  a  trunk  of  Elephas  maximus,  the  Asian  elephant. 
Longitudinal  sheets  of  radial  muscles  radiate  from  the  connective  tissue  around 
the  two  nostrils.  These  muscles  insert  dorsal ly  and  laterally  into  tendinous 
leaves,  not  shown  in  this  figure,  which  in  turn  radiate  and  connect  to  a  long 
tendinous  sheet,  which  is  shown,  covering  the  dorsal  and  lateral  muscle  of  the 
trunk.  The  outer  insertion  of  the  ventral  radial  muscles  is  not  known.  Longi¬ 
tudinal  muscles  run  in  the  channels  formed  dorsal ly  and  laterally  by  the  radi¬ 
ating  muscles.  Some  of  the  longitudinal  muscles  insert  into  the  radiating 
tendons  attached  to  the  radiating  muscles,  and  some  end  in  the  center  of  the 
trunk  where  their  insertion  is  less  obvious. 

Two  layers  of  oblique  muscle,  oriented  in  opposite  directions  and  one 
interior  to  the  other,  run  down  the  ventral  side  of  the  trunk.  Both  of  these 
muscle  masses  are  symmetrical  about  the  sagittal  plane.  The  outer  layer 
inserts  ventral ly  along  the  sides  of  a  long  ventral  tendon  and  runs  laterally 
to  insert  along  the  sides  of  the  dorso-lateral  tendinous  sheet.  The  inner 
layer  inserts  ventrally  into  the  inner  portion  of  the  ventral  tendon,  and  the 
lateral  insertion  of  this  muscle  mass  is  not  known.  The  orientation  of  the 
oblique  muscles  changes  down  the  length  of  the  trunk.  Proximal ly  the  outer 
layer  is  almost  transverse  in  orientation,  but  it  gradually  changes  so  that 
distal ly  the  fibers  run  almost  longitudinally.  The  inner  layer  is  exactly  the 
reverse.  Proxlmally,  the  inner  oblique  muscle  fibers  run  almost  longitudinally 
and  distally  they  run  almost  transversely.  Both  the  dorso-lateral  and  ventral 
tendons  are  separated  proximally  from  the  skin  by  a  layer  of  loose  connective 
tissue;  further  down  the  trunk  they  are  more  intimately  connected  to  the  skin. 
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Figure  1  Cross  section  of  an  elephant's  trunk  showing  the  distribution  of 
muscles  (based  on  information  taken  from  Boas  and  Paul  1 i  (1908) 
by  Brad  Smith) 
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We  took  16  mm  movies  of  two  mature  Asian  elephants  at  the  National  Zoo  in 
Washington,  D.C.  This  involved  filming  the  trunks  while  the  elephants  lifted 
payloads  at  the  tips  of  their  trunks.  These  elephants  had  trunks  about  1.5  m 
long  when  in  a  fully  extended,  relaxed  state.  The  nominal  outside  diameters  of 
these  trunks  were  about  32  cm  at  the  head  and  8  cm  at  the  tip. 

For  our  analysis,  we  used  both  data  from  films  and  data  obtained  from 
drawings  of  cross-sections  of  an  elephant's  trunk  from  the  Boas  and  Paul  1 i 
treatise.  A  composite  of  such  cross-sections  is  shown  in  Figure  1.  Based  on 
the  position  of  radial  muscles  below  the  nostrils,  the  shape  of  the  nostrils, 
and  the  relative  thickness  of  the  oblique  muscles  in  cross-section,  we  deduced 
that  this  particular  transverse  cross-section  was  located  at  position  A-A  along 
the  trunk  as  shown  in  scale  drawings  of  photographs.  Figures  2.  The  length 
scales  shown  in  Figures  1  and  2  are  those  for  the  Asian  elephant  that  we 
filmed.  The  important  assumption  we  used  in  deriving  the  particular  numerical 
results  that  follow  is  that  the  muscle  mass  distribution  for  section  A-A  in 
Figure  2  is  geometrically  similar  to  that  of  Figure  1. 

LIFTING  EXPERIMENTS 

In  preparation  for  a  typical  experiment  involving  a  live  Asian  elephant, 
white  zinc  oxide  spots  were  first  painted  along  the  midlateral  line  on  the  side 
of  the  trunk,  as  shown  in  Figure  2.  Each  white  spot  was  located  along  the 
trunk  at  a  distance  from  the  spot  nearest  the  head.  The  trunk  was  photo¬ 
graphed  in  what  appeared  to  be  Its  longest  and  fully  relaxed  state,  in  a  near 
vertical  position. 

The  sequence  of  events  involving  the  elephant  lifting  a  payload  with  the 
tip  of  Its  trunk  was  recorded  using  a  Canon  Scopic  16  nm  movie  camera  situated 
about  12  m  directly  to  the  side  of  the  animal.  Trunk  motion  was  in  the 


Figure  2  Trunk  contracted  with  payload  on  ground  (a);  payload  at  an 

Intermediate  position  (b);  and  payload  at  Its  highest  position 
(c).  Traced  from  projected  movie  Images. 
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vertical  plane,  parallel  to  the  side  of  the  animal  and  perpendicular  to  the 
plane  of  movie  images.  The  sequence  of  events  was  as  follows. 

1.  The  trunk  tip  looped  around  the  rope  attached  to  a  payload  mass  of 
26.2  kg.  With  the  payload  still  on  the  ground,  the  entire  trunk  shortened  and 
thickened  before  lifting,  resulting  in  the  configuration  of  Figure  2a.  The 
average  longitudinal  contracted  strain  eu  was  deduced  by  comparing  the  loca¬ 
tions  of  the  trunk  spots  in  the  relaxed  and  contracted  states.  That  is,  each 
spot  moved  from  its  relaxed  position  if  by  an  amount  Aif  as  depicted  by  the  in¬ 
sert  of  Figure  3.  The  results  are  plotted  as  circles  in  Figure  3,  from  which 
the  slope  or  initial  strain  before  lifting  was  deduced  by  a  least  squares 
straight  line  fit  as  eu  *  -24.4%,  a  nearly  uniform  value  along  the  trunk. 

2.  The  coiled  trunk  tip  stretched  a  bit  while  the  elephant  lifted  its 
head  and  dragged  the  payload  along  the  ground. 

3.  The  payload  was  lifted  by  an  upward  motion  of  the  head,  accompanied  by 
a  further  shortening  of  the  trunk,  and  a  further  curling  motion  of  the  whole 
trunk.  Figures  2b  and  2c  depict  the  trunk  with  the  payload  lifted  to  an  inter¬ 
mediate  and  the  highest  position.  The  average  value  of  eu  for  these  last  two 
positions  based  on  the  data  of  Figure  3  were  -31.5%  and  -33.0%  respectively. 

A  careful  study  of  enlarged  prints  from  the  film  showed  that  the  radius  of 
curvature  p  traced  by  the  white  dots  in  the  plane  of  motion  changed  signifi¬ 
cantly.  The  values  of  p  and  the  corresponding  uncertainties  in  measurement  ± 

Ap  at  section  A-A  of  the  trunk  were  measured  as  p0  *  1.18  m  ±  0.06  m  and  pj  * 
1.67  m  ±  0.25  m,  corresponding  respectively  to  the  intermediate  and  highest 
positions  of  the  trunk  shown  in  Figure  2.  The  analysis  that  follows  relates 
these  bending  strains  to  the  radius  of  curvature. 
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Figure  3  Straight  lines  Indicate  uniform  trunk  contraction  (constant 

gradient  of  deformation)  for  the  three  conditions  of  Figure  2 
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MATHEMATICAL  MODEL  AND  ANALYSIS 


The  purpose  of  this  section  is  twofold:  to  define  a  rational,  static, 
mathematical  model  of  the  trunk  muscle  system  that  supports  both  the  trunk 
load  and  payload  and  to  calculate  the  apparent  elastic  modulus  of  the  trunk  at 
section  A-A  when  the  trunk  is  between  the  two  positions  shown  in  Figure  2b  and 
2c.  A  few  simplifying  assumptions  are  needed.  First,  the  longitudinal  strain 
ex  due  only  to  trunk  curvature  p  has  a  linear  distribution  from  dorsal  to 
ventral.  That  is 


where  Figure  4  depicts  this  strain  at  position  y  from  the  neutral  axis  of 
bending.  This  linear  strain  distribution  superimposed  on  the  uniform  con¬ 
tractual  strain  eu  from  muscle  pretension  gives  the  total  strain  at  the  trunk 
cross  section.  Second,  in  comparison  to  this  total  strain  distribution,  it  is 
reasonable  to  assume  that  3  r^y  “0,  where  By  is  the  normal  strain  along 
the  y  direction  and  yXy  is  the  transverse  shear  strain  due  to  changes  in  the 
bending  moment  and  curvature  along  the  trunk.  These  assumptions,  which  are 
consistent  with  elementary  beam  theory,  lead  to  an  expression  for  the  normal 
strain  cq  along  the  length  of  each  oblique  muscle  inclined  at  an  acute  angle  0 
with  the  longitudinal  axis  (Figure  3  and  4)  or 

eg  *  excos2e  (2) 

Equation  (2)  follows  from  the  strain  transformation  equations  derived,  for 
instance,  by  Timoshenko  and  Goodier  (1951),  where  ty  and  r^y  are  zero. 

The  third  assumption  is  that  the  longitudinal  tensile  stress  in  each 
muscle  mass  is  linear  with  its  longitudinal  contracted  strain,  or 

°x  ■  -  Ea  ®x*  °e  *  ’  Ea  e0 

where  ox  and  o$  are  the  normal  stresses  for  the  longitudinal  muscles  and 
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oblique  muscles,  respectively.  Here,  Eg  is  the  apparent  elastic  modulus  of 
the  muscle  mass,  which  is  assumed  to  be  the  same  for  each  muscle  type. 

The  fourth  assumption  is  that  only  the  muscles  carry  the  stresses  set  up 
by  the  trunk's  own  mass  and  by  the  payload.  Thus,  the  tissue  surrounding  the 
muscles,  although  following  the  linear  strain  distribution,  has  a  negligible 
stiffness  compared  to  the  stiffness  of  the  activated  muscles.  A  consequence 
of  these  assumptions  is  that  the  neutral  axis  in  bending  does  not  shift  for 
small  changes  in  trunk  curvatures.  The  accuracy  of  these  assumptions  is  not 
known. 


At  a  typical  cross  section  such  as  A-A  in  Figure  4,  the  condition  of 

force  equilibrium  in  the  x  direction,  or  perpendicular  to  this  transverse 

section,  must  be  met.  This  equilibrium  condition  establishes  the  location  of 

the  neutral  axis  y*0,  or  the  plane  along  which  the  muscle  strain  due  to  trunk 

curvature  alone  is  always  zero.  Force  equilibrium  is  expressed  as: 

/  ox  dAx  +  /  oe  cos2e  dA0  *  0  (4) 

Ax  A0 

where  Ax  is  the  total  transverse  cross  sectional  area  for  the  longitudinal 
muscles  and  Aq  Is  the  total  observed  cross  sectional  area  of  the  oblique 
muscles  projected  on  the  same  transverse  cross  section  (the  x*0  plane).  From 
Figure  5,  it  is  observed  that  for  the  oblique  muscles  the  cross  section  area 
element  normal  to  o6  is  dAocos0,  that  the  corresponding  force  is  o8  cose  dA0, 
and  the  component  of  this  force  projected  along  x  is  thus  oq  cos2e  dA0.  which 
leads  to  the  second  integral  on  the  left  in  Equation  (4).  With  Equations  (1), 
(2)  and  (3),  Equation  (4)  is  rewritten  as 


/  ydAx  ♦  /  y  cos*e  dA0  *  0 
Ax  A0 

A  practical  way  to  locate  the  neutral  axis  is  as  follows. 


(5) 

As  shown  in 


Figure  4,  the  neutral  axis  at  y  ■  0  is  located  at  a  distance  z  from  a  datum 
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Figure  4.  Definitions  of  trunk  geometry  and  longitudinal  strain. 
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axis  chosen  as  the  ventral  edge.  If  the  centroid  of  the  area  for  each  muscle 
1  of  area  Ax^  or  A0^  Is  located  at  y  ■  y^,  then  the  centroidal  distance  z-j 
from  the  datum  axis  is  defined  as 


Zi  *  z  -  y1  (6) 
When  Equation  (6)  is  combined  with  Equation  (5),  with  the  differential  areas 
replaced  by  the  areas  Ax^  and  A^ ,  and  the  integrals  replaced  by  finite  sums, 
the  result  is  a  practical  formula  for  locating  the  neutral  axis.  That  is 


w  ^  +  w  2iAoicos4®i 
R  R  7 

I  Axi  +  l  AoiCOS  61 
1-1  X1  1-1  01  1 


where  there  are  N  longitudinal  muscles  and  M  oblique  muscles  at  a  cross 
section.  Care  must  be  taken  to  use  the  correct  angles  for  the  oblique  muscle 
areas.  For  instance,  for  section  A-A  of  Figure  5,  e  -  0,  *  29.1  deg  for  the 
inner  oblique  muscles  and  0  =  *  34.1  deg  for  the  outer  oblique  muscles. 

Having  located  the  neutral  axis,  the  bending  strain  distribution  can  be 
calculated  from  Equation  (1)  If  the  radius  of  curvature  is  known.  Of  interest 
here  is  the  change  in  the  bending  strain  distribution  between  two  loading 
states  at  section  A-A.  For  instance,  the  difference  in  longitudinal 


muscle  strain  Aex  for  a  trunk  configuration  in  its  highest  position  (Figure 
2c)  for  p  *  pi  and  in  its  intermediate  position  (Figure  2b)  for  p  *  p0  is 

Ac  *  -  X—  -  ( -  ^— )  *  -  ({ 


where 


po  -  P1 


At  a  fixed  location  y,  the  change  in  strain  Acq  in  a  typical  oblique  muscle 
will  always  be  less  than  Aex,  or  from  Equation  (2) 


Accompanying  these  strain  changes  Aex  and  Ae0  during  lifting  are 
corresponding  changes  in  the  normal  stress  distributions,  Aox  and  Ao0.  These 
are  the  muscle  reactions  in  response  to  shifts  in  two  types  of  external  loads 
causing  moments  at  section  A-A.  One  of  these  loads  is  due  to  the  trunk  mass 
segment  mt  (below  section  A-A  to  the  tip)  which  causes  a  differential  moment 
of  mtg(di-d0).  The  other  load  is  due  to  the  payload  mass  nip  which  causes  a 
differential  moment  mpg(ei-e0).  The  load  offset  distances  d0,  dj.,  e0  and 
are  defined  in  the  free  body  sketches  of  Figure  6.  For  moment  equilibrium  at 
Section  A-A,  it  follows  that 

mtg(d1-d0)  +  mpgUi-eJ  =  /  y  AoxdAx  +  /  y  Ao0  cos2edA0  (11) 

Ax  Ao 

Here  Aox  dAx  and  Ao0  cos29  dA0  are  the  changes  of  the  x-directed  forces  in  the 
longitudinal  and  oblique  muscles,  respectively,  which,  when  multiplied  by 
their  moment  arms  y  and  integrated  over  section  A-A,  leads  to  Equation  (11). 
Invoking  the  linear  relationships  of  Equations  (3),  along  with  Equations  (2), 
(8),  and  (10),  the  differential  stresses  become 


l0x  ■  -  Ea  “x  ■  Ea  J 


4oe  ■  Ea  “a  *  Ea  -  cos  6 

P 


(12) 

(13) 


With  Equations  (12)  and  (13),  Equation  (11)  can  be  solved  for  the  apparent 
modulus,  or 


Ea  -  7  [V(dl'V  +  U4) 

where  p  is  given  by  Equation  (9)  and  I,  the  second  moment  of  the  muscle  area 
with  respect  to  the  neutral  axis  (y*0),  is 


With  a  knowledge  of  the  muscle  morphology  as  shown  in  Figure  1,  along  with 
trunk  curvatures  and  the  loading  terms,  the  apparent  modulus  Ea  may  be 
calculated  from  Equations  ( 14)  and  (15).  Then  the  differential  stress  and 
strain  distributions  due  only  to  a  curvature  change  may  be  calculated  from 
Equations  (12)  and  (13).  The  appropriate  measurements  and  typical  numerical 
results  are  now  discussed. 

NUMERICAL  RESULTS  AND  DISCUSSION 

To  obtain  the  desired  numerical  results,  a  clear  and  properly  scaled  out¬ 
line  of  the  longitudinal  and  oblique  muscle  areas  at  Section  A-A  is  needed. 
Such  a  drawing,  deduced  from  the  Boas  and  Paulli  treatise  (1908)  is  shown  in 
Figure  7.  This  muscle  distribution  and  the  measured  data  listed  in  Table  1 
were  used  to  calculate  the  apparent  trunk  modulus.  The  geometry  of  Figure  7 
was  scaled  to  match  the  cross  section  A-A  of  our  subject  animal's  trunk,  for 
which  the  section  height  D  (Figure  4)  was  measured  from  photographs  as  24  cm. 
Of  the  150  muscle  areas  shown  in  Figure  7,  78  are  longitudinal,  51  are  inner 
oblique  (8  *  29.1  deg),  and  21  are  outer  oblique  (8  »  34.1  deg).  Using  a  fine 
grid  overlay,  each  muscle  area  or  projected  area  and  its  respective  centroidal 
distance  to  the  datum  plane  was  calculated.  The  total  areas  Ax  and  A0  are 
given  in  Table  2.  With  the  values  of  8  measured  from  drawings  in  the  Boas  and 
Paulli  treatise  (1908),  z  was  calculated  from  Equation  (7)  as  12.0  cm.  Since 
D  *  24  cm,  the  neutral  axis  lies  at  the  mid-height  of  this  cross  section.  The 
value  of  I  was  then  calculated  by  numerically  integrating  Equation  (15), 
yielding  I  *  4270  cm*.  Of  course,  these  numerical  results,  although  precisely 
calculated  to  three  significant  figures,  are  only  as  accurate  as  the  artist’s 
original  drawings  of  the  muscle  distributions. 
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Idealized  muscle  distribution  used  for  calculations 
(adapted  from  Boas  and  Paulli,  1908) . 


D  ■  24  cm 


d0  s  5.0  cm 
di  =  8.0  cm 
e0  =  9.29  cm 
ej  =  13.27  cm 

Ai »  Aii 
mp  =  26.2  kg 

mt  a  16.37  kg 
M  *  78 
N  =  72 

Y  =  1180  kg/m^ 
e  *  29.1  deg 
0  =  34.1  deg 

p0  - 

PI  =  167 
Ap0  =  i  6  cm 
API  *  i  25  on 


trunk  height  (y  direction)  at  section  A-A 
offset  distance  of  (intermediate).  Fig.  6a 

offset  distance  of  (highest).  Fig.  6b 

offset  distance  of  mp  (intermediate).  Fig.  6a 

offset  distance  of  mp  (highest).  Fig.  6b 

spot  locations,  plotted  in  Fig.  3 
mass  of  payload 

mass  of  trunk  from  tip  to  section  A-A 
number  of  longitudinal  muscles  at  section  A-A 
number  of  oblique  muscles  at  section  A-A 
mass  density  of  trunk 

inclination  of  inner  oblique  muscles.  Fig.  5 

inclination  of  outer  oblique  muscles.  Fig.  5 

radius  of  curvature  of  section  A-A,  Fig.  6A 
radius  of  curvature  at  section  A-A,  Fig.  6b 
uncertainty  of  p0  measurement 

uncertainty  of  pi  measurement 
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Table  2  -  Summary  of  Calculated  Data 


Symbol 
and  Value 


Meanirn 


Equation 
or  Figure 


Ax  *  54.6  cm^ 

area  of  all  longitudinal  muscles  on 
x=0  plane 

Fig.  7 

A0  *  185.7  cm2 

projected  area  of  all  oblique  muscles  on 
xa0  plane  (inner  oblique:  160  cm2;  outer 
oblique:  25.7  cm2) 

Fig.  7 

Ea  *  1.33x11)6  N/m2 

apparent  modulus  of  trunk 

Eq.  (14) 

AEa  *  ±  0.532x106  N/m2 

uncertainty  in  Ea 

Eq.  (16) 

I  *  4270  cm4 

area  moment  of  all  muscles 

Eq.  (15) 

z  *  12.0  cm 

location  of  neutral  axis 

Fig.  4 

aox  *  3.22xl04  N/m2 

max.  stress  change,  longitudinal  muscle 

Eq.  (12) 

Aog  «  2.23xl04  N/m2 

max.  stress  change,  inner  oblique  muscle 

Eq.  v13) 

Aoq  =  2.37xl04  N/m2 

max.  stress  change,  outer  oblique  muscle 

Eq.  (13) 

Aex  *  0.0242 

max.  strain  change,  longitudinal  muscle 

Eq.  (8) 

Aeg  *  0.0168 

max.  strain  change.  Inner  oblique  muscle 

m 

• 

M 

o 

Acg  *  0.0178 

max.  strain  change,  outer  oblique  muscle 

Eq.  (10) 

eu  »  -0.244 

uniform  strain,  lowest  position 

Fig.  3 

eu  *  -0.315 

uniform  strain,  intermediate  position 

Fig.  3 

eu  *  -0.330 

uniform  strain,  highest  position 

Fig.  3 

p  *  -402  cm 

differential  curvature 

Eq.  (9) 

The  external  moments  at  section  A-A  due  to  the  payload  and  the  trunk  mass 
were  measured  for  the  trunk  in  two  different  positions:  the  intermediate  and 
the  highest  positions  shown  in  Figures  2  and  6.  The  reason  for  choosing  these 
two  positions  was  that  the  uniform  component  of  longitudinal  trunk  strain  was 
nearly  the  same  for  both  positions  (eu  was  0.315  and  0.330,  respectively). 
Thus,  the  radii  of  curvature  p0  and  pj  as  measured  from  photographs  did  not 
need  to  be  corrected  to  account  for  changes  in  eu.  The  importance  of  highly 
accurate  measures  of  these  radii  to  an  accurate  calculation  of  the  apparent 
modulus  will  be  discussed  below. 

The  quantities  defining  the  differential  moment  between  the  selected 
trunk  positions  are  enclosed  in  the  square  brackets  of  Equation  (14). 
Measurements  were  made  for  each  of  these  quantities.  The  payload  mass  mp  was 
26.2  kg.  The  offset  distances  e0  and  e^  for  this  load,  as  well  as  the  volume 
segment  of  the  trunk  from  tip  to  section  A-A,  were  deduced  from  photographs 
(see  Figure  6).  The  centroidal  distances  dj  and  d£  for  this  volume  segment 
were  calculated  assuming  that  the  volume  was  a  series  of  right,  truncated 
cones.  The  volumes  of  the  two  nostrils,  deduced  from  the  Boas  and  Paulli 
treatise,  were  treated  as  voids.  The  product  of  this  volume  segment  and  the 
trunk's  mass  density  then  gave  the  value  m*  *  16.37  kg.  Our  measurements  of 
the  mass  density,  made  on  selected  samples  from  different  areas  of  a  dead 
elephant's  trunk,  yielded  an  average  mass  density  of  1180  kg/m-*,  with  a 
maximum  variation  of  8  percent  among  samples. 

The  remaining  quantity  needed  to  calculate  Ea  given  by  Equation  (14)  is 
the  differential  curvature  p.  This  is  defined  by  Equation  (9)  in  terms  of  the 
radii  of  curvatures  p0  and  p^  at  section  A-A.  Accurate  measures  of  these 
radii  were  more  difficult  to  deduce  from  enlarged  photographs  than  the  data 
previously  discussed.  More  closely  spaced  and  more  carefully  aligned  white 


spots  along  the  trunk  would  lead  to  better  measurements.  In  the  present  case, 
P0  =  118  ±  6  cm  and  p\  *  167  ±  25  cm,  showing  an  uncertainty  of  5  percent  and 
15  percent,  respectively.  For  the  nominal  values  of  the  radii.  Equation  (9) 
gave  p  *  -402  cm. 

Using  the  measured  data  discussed  above  and  summarized  in  Table  1,  the 
apparent  modulus  was  calculated  from  Equation  (14),  or  Ea  =  1.33  x  10®  N/m2. 

It  is  interesting  to  note  that  this  modulus  compares  favorably  with  the 
tangent  modulus  of  human  Sartorius  muscle  (1.77  x  10®  N/m2)  and  dog  Rectus 
abdominis  muscle  (4.0  x  10®  N/m2)  at  prestrains  of  60  to  65  percent,  as 
estimated  from  data  presented  by  Yamada  and  Evans  (1970).  Also,  this  modulus 
agrees  within  10  percent  with  Young's  modulus  of  a  commercially  available, 
synthetic  rubber.  Type  A6UB*. 

Further  calculations  lead  to  the  differential  strain  distributions  in  the 
longitudinal  muscles,  Aex,  and  in  the  oblique  muscles,  Aeq,  and  their 
corresponding  differential  stresses,  Aax  and  Aoe.  The  maximum  value  for  each 
of  these  quantities  was  calculated  from  Equations  (8),  (10),  (12)  and  (13) 
based  on  the  maximum  distance  y  extending  from  the  neutral  axis  to  the 
centroid  of  the  most  outlying  muscle.  The  numerical  results  summarized  in 
Table  2  reveal  that  the  maximum  differential  muscle  strains  range  from  about 
1.7  to  2.5  percent,  and  the  maximum  differential  stresses  range  from  about  2.2 
x  1U4  to  3.2  x  104  N/m2.  In  fact,  at  any  chosen  distance  y  from  the  neutral 
axis,  calculations  showed  rather  narrow  ranges  for  the  differential  strains 
and  differential  stresses  among  the  three  types  of  muscles.  It  is  apparent 
that  the  muscles  do  a  good  job  in  sharing  the  load  at  this  cross  section. 

Given  the  validity  of  the  mathematical  model,  for  which  the  main 
assumptions  were  the  linearity  of  strain  over  the  cross  section  and  the  linear 
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stress-strain  law,  the  accuracy  of  the  numerical  results  were  found  to  depend 
strongly  on  the  accuracy  of  p0  and  pj.  These  were  the  most  uncertain  of  all 
of  the  measured  parameters.  The  analysis  of  Kline  and  McClintock  (1953)  may 
be  used  to  show  just  how  sensitive  the  calculation  of  Ea  is  to  the  measured 
uncertainties  Ap0  and  Apj  of  p0  and  pj,  respectively.  The  corresponding 
uncertainty  in  Ea,  or  AEa  for  these  two  parameters  only  is  given  by 

>E_  \  2  /3Ea  \211/2 

iEa  '  lpo)  +  4pl)  j  (16> 

When  p  from  Equation  (9)  is  substituted  into  Equation  (14)  and  the  result  is 
differentiated  according  to  Equation  (16),  the  result  is 


Vpi' 


4po2  +  po4 


<) 


2\  1/2 


With  the  nominal  values  of  the  radii  and  their  deviations  as  listed  in  Table 
1,  the  numerical  value  of  AEa/Ea  calculated  from  Equation  (17)  is  0.40.  Thus, 
for  relatively  small  uncertainties  In  radii  (5  percent  and  15  percent),  the 
uncertainty  in  the  apparent  modulus  Is  40  percent,  where  Ea  *  1.33  x  106  H/m2 
±  0.532  x  106  H/m2. 

Further  studies  are  now  underway  that  involve  microscopic  studies  of 
muscle  morphology  at  different  cross  sections  in  a  dead  elephant's  trunk,  and 
also  the  direct  measurements  of  tensile  properties  of  connective  tissue.  As 
further  data  becomes  available,  the  present  mathematical  model  may  be  refined. 
Thus,  more  exact  connections  may  be  deduced  among  the  trunk's  system 
parameters  Including  its  external  loading,  its  overall  shape  during  lifting, 
its  muscle  distribution,  and  the  mechanical  properties  of  its  muscles.  The 
results  may  provide  clues  to  future  mechanical  designs  of  continuously 
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NOMENCLATURE 


A0  total  area  of  oblique  muscles  projected  on  the  x=0  plane 

A0i  projected  area  of  1-th  oblique  muscle  on  the  x=Q  plane 

Ax  total  cross  sectional  area  of  longitudinal  muscles  in  the  x=0  plane 

Ax-j  cross  sectional  area  of  i-th  longitudinal  muscle  on  the  x=0  plane 

D  trunk  height  in  y  direction  at  section  A-A 

d0,d^  moment  arms  of  m^,  Figure  6 

Ea  apparent  modulus  of  trunk 

A£a  uncertainty  in  Ea 

e0,e-j  moment  arms  for  mp.  Figure  6 

g  gravitational  constant 

I  second  moment  of  muscle  cross  sectional  areas,  Equation  (15) 

t-j  distance  of  the  i-th  spot.  Figure  3 

M  number  of  oblique  muscles  at  a  cross  section 

mp  mass  of  pay  load  lifted  at  trunk  tip 

mt  mass  of  trunk  from  tip  to  section  A-A,  Figure  6 

N  number  of  longitudinal  muscles  at  a  cross  section 

x  longitudinal  coordinate  along  trunk's  neutral  axis 

y,y^  transverse  coordinates  measured  from  neutral  axis.  Figure  4 

Zi  transverse  coordinate  of  1-th  muscle.  Figure  4 

z  location  of  neutral  axis  from  ventral,  Figure  4 

YXy  transverse  shear  strain 

ex  normal  strain  in  longitudinal  muscles 

Cy  normal  strain  along  y 

ce  normal  contraction  strain  in  oblique  muscles 

eu  uniform  trunk  contraction  strain 

a£.i  displacement  of  1-th  spot  along  longitudinal  axis 


Aex  differential  longitudinal  strain  due  to  trunk  curvature  change 

Ae0  differential  strain  in  oblique  muscles  due  to  trunk  curvature  change 

e  acute  inclination  angle  of  oblique  muscles  with  the  x-axis 

Q-j  value  of  6  for  the  i-th  oblique  muscle 

p  radius  of  curvature  of  trunk's  neutral  axis 

p0  value  of  p  at  section  A-A,  intermediate  load  position 

Ap0  uncertainty  in  p0 

Pl  value  of  p  at  section  A-A,  highest  load  position 
Apj  uncertainty  in  pj 

p  differential  radius  of  curvature.  Equation  (9) 

ox  normal  tensile  stress  in  longitudinal  muscles 
Oy  normal  stress  in  transverse  direction 
00  normal  tensile  stress  in  oblique  muscles 

aox  differential  stress  corresponding  to  Aex 

ao6  differential  stress  corresponding  to  Ae0 
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II.  MECHANICS  OF  A  CONTINUOUS  MANIPULATOR  MADE  OF  A 


NONLINEAR,  COMPOSITE  MATERIAL 


ABSTRACT 


This  research  deals  with  the  analysis  of  light  weight,  flexible  and 
fast  moving  robotic  arms.  A  typical  arm  consisting  of  hollow  cylindrical  ele¬ 
ments,  bends  and  coils  around  a  load  and  moves  it  to  the  desired  position  by 
varying  the  internal  pressures  in  the  individual  elements.  The  theoretical 
derivations  relating  arm  positions  and  loading  incorporate  nonlinear  material 
properties  and  large  deflections.  A  learning  program  is  developed  that  stores 
information  about  arm  motion,  information  that  is  used  in  subsequent  manipula¬ 
tions  to  achieve  efficient  lifting  scenarios.  Typical  numerical  results  for 
arm  designs  will  be  verified  experimentally  in  the  near  future. 


INTRODUCTION 


Robots  are  programmed  machines  designed  to  carry  out  desired  tasks  for 
humans  [1,2].  Their  use  in  the  industry  is  rapidly  increasing  and  extensive 
research  is  underway  with  the  aim  of  improving  their  capabilities.  Robots  can 
be  used  to  efficiently  carry  out  routine  tasks  without  making  the  mistakes 
that  a  human  would  because  of  fatigue.  Also  robots  can  be  used  in 
environments  which  are  harmful  to  the  health  of  humans. 

Mecnanically  most  of  the  present  day  robots  have  arms  consisting  of 
rigid  members  connected  by  hinged  or  pivoted  joints.  To  move  such  an  arm 
precisely,  the  movements  of  its  various  members  have  to  be  coordinated.  Such 
robots  are  slow-acting  because  of  the  inertia  of  the  members  and  the  time 
required  to  compute  and  control  position. 

A  different  approach  using  a  light-weight,  highly  flexible,  fast 
moving  arm  is  made  in  this  research.  Such  an  arm  is  analogous  to  an 
elephant's  trunk  which  has  been  analyzed  recently  [3],  It  was  observed  that 
trained  elephants  can  efficiently  coil  their  trunks  around  heavy  loads  and 
lift  them  to  desired  positions.  The  present  study  simulates  the  action  of  an 
elephant's  trunk  within  the  limits  of  structural  feasibility.  The  robotic  arm 
is  a  continuous  manipulator  of  rubber  or  polymeric  material,  consisting  of 
hollow  cylindrical  elements  joined  end  to  end.  The  elements  are  selectively 
reinforced  so  that  they  bend,  twist  or  extend  when  pressurized.  Theoretical 
derivations  and  computer  programs  are  written  to  evaluate  the  internal 
pressure  required  in  each  of  the  elements  such  that  the  arm  can  coil  around  a 
load  and  lift  it  to  a  desired  position. 

This  research  also  deals  with  another  important  aspect  of  robotics. 

As  machines  cannot  reason  or  think  like  humans,  the  robot  has  to  be  programmed 
to  anticipate  the  various  options  or  remedial  measures,  in  advance  of  motion. 
Thus  some  aspects  of  Artificial  Intelligence  are  Incorporated  into  this 
research.  A  new  Learning  Program  evaluates  the  end  moments  to  be  applied  to 


each  of  the  elements  in  order  to  lift  and  position  the  end  load.  This 
program's  database  is  updated  appropriately  each  time  the  program  is  run, 
thereby  increasing  the  apparent  knowledge  and  efficiency  of  the  program  in 
performing  load  positioning  calculations. 

The  theoretical  derivations  for  the  mechanical  behavior  of  the 
pressurized  elements,  the  programming  logic,  and  some  numerical  results  for  a 
typical,  flexible  robotic  arm  are  discussed. 

ANALYSIS  OF  A  TYPICAL  ELEMENT 

The  geometry  of  a  typical  element  is  shown  in  Figure  1.  It  is  a  thin 
walled  cylinder  with  the  ends  sealed  and  with  reinforcement  on  the  underside. 
The  element  is  modelled  as  a  cantilever  beam  of  a  nonlinear  material.  As  in 
[4],  finite  deflections  are  considered.  Given  the  applied  end  moment  due  to 
internal  pressure,  the  horizontal  and  vertical  deflections  of  the  end  of  the 
beam  are  evaluated.  The  radius  of  curvature  of  the  deflected  beam  is  also 
computed.  In  this  analysis,  [5]  was  also  used. 

The  hollow,  cylindrical  section  is  made  out  of  a  synthetic,  nonlinear 
material  for  which  the  stress  strain  relationship  is  as  follows 

a  =  Ae  +  Bel/n 

where  A,  B  and  n  are  material  constants. 

The  synthetic  material  can  withstand  high  strains  (up  to  200%)  without 
failure.  The  reinforcement  is  made  of  high  strength  steel  shim  or  synthetic 
fibers  with  high  tensile  strength  (Kevlar,  for  instance)  as  shown  in  Figure  2. 
The  reinforcement  provides  the  strength  to  the  structure.  The  neutral  axis 
for  the  composite  element  passes  through  the  reinforcement  because  its  stiff¬ 
ness  is  much  greater  in  magnitude  than  that  of  the  synthetic  cylindrical 
section.  For  this  reason  with  the  application  of  internal  pressure  and  thus 
an  end  moment,  the  element  bends  as  shown  in  Figure  1.  The  reinforcement  is 
assumed  to  be  a  linear  material. 
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pressurization. 


Because  the  cylindrical  part  of  the  composite  element  is  made  up  of  a 
nonlinear  composite  material,  the  position  of  the  neutral  axis  varies  and  has 


to  be  computed  for  each  stress  distribution  across  the  section.  Given  an 
applied  end  moment  Ma,  the  following  iterative  steps  are  required  to  calculate 
the  position  of  the  neutral  axis. 

a.  Choose  an  initial  position  of  the  neutral  axis. 

b.  Compute  the  radius  of  curvature  p  for  the  deflected  element  for 
the  chosen  position  of  the  neutral  axis. 

c.  Compute  the  moment  of  resistance  Mc  of  the  section  for  the 
computed  p  and  the  chosen  position  of  the  neutral  axis. 

d.  If  the  following  inequality  holds  for  a  specified  tolerance  limit 
e 


then  the  iterations  stop.  Otherwise  a  new  position  of  neutral 
axis  is  chosen  and  steps  (b)  to  (d)  are  repeated.  See  Appendix  A 
for  the  details  of  calculation. 


Derivation  of  End  Deflection  Formulae 


After  the  position  of  the  neutral  axis  and  the  radius  of  curvature 
have  been  determined  for  the  given  end  moment,  the  end  deflections  are 
computed  in  the  following  manner. 

We  know  the  moment  curvature  relationship  as  given  by  equation  (A. 21 
in  Appendix  A).  The  equation  is  repeated  here  as  follows 


Mc  •  (±>«  *  <i)6  *  (i)1/n 


On  rearranging  the  terms  we  yet 


p 


U  +  Cy  (x)]V^ 


Equatiny  Equations  (1)  and  (2)  we  get 

_ y"(x)  _  .  Mc  - 

a  +  [y'(x)]2}3/2  a+e 

Equation  (3)  is  a  second  order  differential  equation  which  can  be  integrated 
directly  with  respect  to  x  to  yield 


(Mc  -  5) 


_ y  (x )  _ 

U  +  [y  (x)]2)17^  (a'+  3)  *  +  Cl 


The  constant  of  integration,  c\,  can  be  evaluated  using  the  fact  that  y  (0)  = 
0.  Thus  ci  =  0  and 


jlM. 


(Mc  -  O 


a  +  Cy’(x)]2}1/2  (a  +  e) 


Solving  for  y  (x)  yields 


,  (M  -  Ox 

y  (*)  »  C  (4) 

7(a  +  S)Z  -  (Mc  -  02x2 
By  integrating  once  more  one  obtains 

(at  6)  -  [(a  t  6)2  -  (M  -  S)2  x2]1/2 

*(x> - -  (5) 

where  the  constant  of  integration  has  been  evaluated  using  the  condition  that 
y(0)  =  0. 


To  evaluate  the  horizontal  deflection,  6^ ,  at  the  free  end  of  the 
beam,  as  shown  in  Figure  3,  the  equation  for  the  arc  length  is  used.  The 
total  length  then  is 


l  *  L  +  Cy'(x)]2 


dx 


From  Equations  (5)  and  (6)  we  obtain 

L-6. 


L  -  L 


a  +  S 


7(a  +  e)2  -  (Hc  -  5)2x2 


dx 


This  integral  can  be  evaluated  using  the  trignometric  substitution. 

a  +  3 


x  = 


MC-  5 


sin  0 


a  +  8 

dx  =  z - -  cos  0  d  0 


MC  -  5 


From  Equation  (8) 


6  =  s1n'1  I*  -prnsf) 


(H,  -  5) 


when 


x  =  L-6h 
-1 


J(at  x-L-y  -  6  =  s1n  [(L  ’  V  (.  ♦  e)  1 


(11 


Using  Equations  (9)  and  (10)  from  Equation  (7)  we  get 


>  T* 
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ANALYSIS  Of  MOTION  AND  LOAD-CARRYING  CAPACITY  OF  A 
CONTINUOUS  MANIPULATOR  ARM 

The  manipulator  arm  is  made  up  by  joining  a  number  of  the  ‘typical* 
elements,  end  to  end.  The  desired  capability  of  the  arm  determines  the  size 
and  number  of  elements  in  the  arm.  The  last  few  elements,  (3  in  our  model) 
coil  around  the  load,  whereas  the  remaining  elements  help  in  the  lifting  of 
the  load.  Refer  to  Figure  4b. 

Given  an  applied  end  moment,  as  shown  in  the  previous  chapter,  we  can 
evaluate  the  end  deflections,  and  thus  the  end  coordinates  of  the  element  with 
respect  to  the  local  coordinate  axis  system  of  the  element.  These  end 
coordinates  are  transferred  from  local  to  global  coordinate  axis  system  by 
using  a  transformation  matrix. 

Similarly  we  can  carry  out  the  transformations  for  each  of  the 
elements  and  get  the  overall  orientation  of  the  arm  in  the  global  coordinate 
axis  system. 

When  the  end  moments  are  changed,  we  can  by  the  same  procedure, 
compute  the  new  orientation  of  the  arm.  The  self  weight  of  the  elements  is 
neglected  throughout  in  our  derivations. 

The  horizontal  deflection,  at  the  free  end  can  be  obtained  from  Equation 
(11)  as  follows 

„  .  «  (H  -  0 

'h  * L  -  s1n  [L  <12> 

The  vertical  deflection  at  the  free  end  can  be  obtained  by  substituting  (L-6^) 
for  x  in  Equation  (S)  to  get 

6v  =  Smc  -H)  ^  "  cos{L  ( qC+  6) 

Thus  knowing  the  end  moment  values  the  end  deflections  6^  and  <5V  can 
be  computed  from  equations  (12)  and  (13)  respectively,  after  the  position  of 
neutral  axis  and  the  radius  of  curvature  have  been  computed. 
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Derivation  to  Transfer  Local  End  Coordinates  to 
Global  End  Coordinates 


As  the  bending  moment  is  constant  for  an  element,  the  radius  of 
curvature  p  is  constant.  It  follows  from  the  sketch  above  that 


As  shown  in  Figure  5a  for  the  first  element  the  global  and  local  axis 
coincide. 


(14) 


For  the  remaining  elements  the  relationship  between  the  local  and 
global  axis  system  is  as  follows. 
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FIRST  ELEMENT 


V 


Figure  5. 


First  element  showing  coincidence  of  local  and  global 
and  second  element  showing  displaced  local  axis  (b). 
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^0  is  the  angle  which  the  x'  axis,  of  the  local  coordinate  axis  system 
of  the  element,  makes  with  the  x  axis  of  the  global  coordinate  axis 
system. 

X2,y^  are  the  coordinates  of  the  origin  of  the  local  coordinate  axis 
system  in  terms  of  the  global  coordinate  axis  system. 
xl,yi  are  the  coordinates  of  any  point  P  in  the  local  coordinate  axis 
system. 

x,y  are  the  coordinates  of  the  same  point  P  in  the  global  coordinate 
axis  system. 

In  our  case  the  point  P  is  the  end  point  of  the  deflected  element. 

X1  ■  L  -  «h 


yi  =  -  *v 


where 


*  horizontal  deflection 
<SV  *  vertical  deflection 


For  the  next  element 


*  a  f  +  — 

on  oc 


The  subscripts  n  and  c  mean  next  and  current,  respectively.  The  oriyin  for 
the  local  coordinate  axis  system  for  the  next  element,  is  the  end  coordinates 


of  the  current  element 

x2n  *  xc 

(17) 

yzn  *  yc 

Now  using  the  transformation  matrix,  the  end  coordinates  of  the  next 
element  can  be  converted  from  the  local  coordinate  axis  system  to  the  global 
coordinate  axis  system.  Similarly,  this  can  be  done  for  all  the  remaining 
elements. 
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Evaluation  of  End  Moments 

It  is  required  to  compute  the  end  moments  to  be  applied  on  each  of  the 
elements  such  as  to  move  the  load  from  its  initial  position  to  the  final 
desired  position. 

An  example  of  a  manipulator  arm  is  shown  in  Figures  4a  and  4b.  The 
first  five  elements  have  the  shim  (or  reinforcement)  on  the  right  side  whereas 
the  remaining  five  elements  have  it  on  the  left  side.  When  these  elements  are 
pressurized  there  will  be  a  reversal  of  curvature  after  the  fifth  element. 

Also  by  positioning  the  elements  in  this  order,  the  arm  has  added  capability 
in  moving  the  load  over  a  wider  region. 

Initially  the  load  is  assumed  to  have  no  mass.  Then  the  end  moments 
Mp,  to  elastically  bend  the  arm  and  move  the  load  to  its  desired  position  are 
computed. 

Secondly,  after  the  arm  has  reached  its  final  configuration,  we  know 
the  coordinates  of  the  end  point  of  each  element  and  also  the  coordinates  of 
the  center  of  gravity  of  the  load.  Thus  the  additional  moment  due  to  the  load 
Mz,  (the  load  is  now  supposed  to  have  mass)  can  now  be  easily  evaluated. 

The  total  moment  Mt  *  mp  + 

Evaluation  of  Mp 

It  is  observed  that  there  is  no  unique  set  of  Mp,  which  can  shift  the 
load  from  its  initial  to  the  final  position.  The  moments  in  some  elements 
could  be  increased  whereas  the  moments  in  some  other  elements  could  be 
decreased  and  the  load  could  possibly  be  still  in  the  same  position.  To 
overcome  this  problem,  an  iterative  procedure  has  to  be  used.  The  iteration 
steps  are  as  follows. 


a.  The  end  moments  Mp  to  coil  the  last  three  elements  around  the  load 
are  first  evaluated. 

b.  For  the  remaining  elements  initial  guess  moment  Mp  values  have  to 
be  provided.  (In  the  next  chapter  a  better  approach  is  explained). 

c.  Knowing  the  end  moments  Mp  in  each  of  the  elements,  we  can  compute 
the  end  deflections  (6^  and  6V)  and  radius  of  curvature  p,  by  the  method 
discussed  in  Chapter  II.  The  transformation  matrices  can  then  be  used  to  get 
the  overall  orientation  of  the  deflected  arm  (method  is  discussed  earlier  in 
this  chapter). 

d.  The  present  position  of  the  center  of  gravity  of  the  load  is  now 
calculated  as  shown  in  Figure  6b: 

xp  =  Xd  -  (R i  +  Ri  +  H)  cos  (¥) 
yp  3  yd  -  (R*  +  Ri  +  H)  sin  (?) 

Now  the  following  inequalities  must  be  simultaneously  satisfied  for  the 
specified  tolerance  limits  and  respectively. 


|xp  -  Xfl  < 

(18) 

|yp  -  yf 1  <  e2 

(19) 

where  xf  and  yf  are  final  desired  position  coordinates  of  the  center  of 
gravity  of  the  load.  If  the  inequalities  (18)  and  (19)  are  not  satisfied,  the 
moments  in  the  elements  have  to  be  varied  in  such  a  way  that  the  load  moves 
closer  to  the  final  desired  position. 

The  present  position  of  the  load  can  be  in  any  one  of  the  four  regions 
shown  in  Figure  7.  The  change  In  moments  can  be  done  logically  as  follows. 
Refer  to  Figures  6a  and  7. 

1.  When  load  is  to  the  right  and  higher: 

XP  >  Xf 

yp  >  yf 
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*i  .yj  > :  INITIAL  LOAD  POSITION 


Figure  7.  The  four  regions  in  which  the  load  could  be 
currently  positioned. 
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Then  increase  the  moments  in  the  elements  of  the  type  L. 

ii.  When  the  load  is  to  the  left  and  higher 

xp  <  xf 

yp  >  yf 

Then  decrease  the  moments  in  the  elements  of  type  L  and  of  type 

U. 

iii.  When  the  load  is  to  the  right  and  lower 

Xp  >  Xf 

yp  <  yf 

Then  increase  the  moments  in  the  elements  of  type  U  and  of  type 

L. 

iv.  When  the  load  is  to  the  left  and  lower 

xp  <  Xf 

yp  <  yf 

Then  decrease  the  moments  in  the  elements  of  type  L. 

After  making  the  requisite  moment  changes  in  the  elements,  which 
depends  on  the  current  orientation  of  the  arm,  we  continue  the  iterative 
procedure  from  step  (b)  onwards,  till  the  condition  in  step  (d)  is  satisfied. 

Step  (a)  is  carried  out  only  once  as  the  moments  required  to  coil  the 
elements  around  the  load  cannot  change  since  the  dimensions  of  the  load  are 
fixed. 

Evaluation  of 

At  the  end  of  the  iteration  procedure,  for  the  evaluation  of  Mp,  the 
load  will  be  somewhere  near  the  final  desired  position.  In  this  orientation 
of  the  arm,  we  know  the  center  of  gravity  of  the  load  and  the  end  coordinates 
of  each  element.  Knowing  this,  we  can  evaluate  the  additional  moment  due  to 


Since  we  do  not  precisely  know  the  point  of  contact  between  the  load 
and  the  last  3  elements  coiled  around  it,  we  assume  the  worst  case  and  apply 
an  additional  moment  as  given  below. 

Mu  =  Wi  (Ri  +  R*  +  H) 

(For  the  last  3  elements  M$,  is  taken  as  same  and  is  on  the  higher  side) 

We  can  compute  the  additional  moments  in  the  other  elements  from  the 
free  body  diagram  shown  in  Figure  8. 

Mi  =  Wi(Ri  +  Rji  +  H)  cos  (6) 

Mjt  for  element  type  U,  say  element  5  is 

M £  =  Wi(xc  -  xd)  -  Mi 
M*  of  element  type  L,  say  element  6  is 

M*  *  Mi  -  W!(xc  -  xe) 

Similarly  we  can  evaluate  M*,  for  all  the  elements. 

Now  for  each  of  the  elements  we  can  evaluate 

Mt  =  Mp  +  M* 

After  the  end  moments  required  to  be  applied  to  each  of  the  elements 
have  been  evaluated,  the  internal  pressures  required  in  each  of  the  elements 
can  easily  be  computed.  As  shown  in  Figure  9  (an  arm  with  3  elements  only), 
there  are  three  end  moments.  Mi,  Mg  and  M3  and  three  different  pressures 

Pi,  p2  and  p3 
where 

2 

M3  =  P3*R23d3 

2 

M2  *  (p2  -  P3)*R22d2 
2 

Ml  =  (Pi  -  P2)  irR2idi 

where  R21,  R22»  R23  are  the  internal  radii  and  di,  dg,  4$  are  the  distances 
between  the  center  of  pressure  and  the  neutral  axis  for  the  three  elements. 


Relationship  between 


internal  pressures  and  end  moments. 


IMPLEMENTATION  OF  THE  LEARNING  CONCEPT 
A  FORTRAN  program  has  been  written  to  evaluate  the  end  moments 
required  in  each  of  the  elements,  to  move  the  load  to  its  desired  position. 

As  mentioned  in  the  previous  chapter,  apart  from  other  input,  initial  guess 
end  moment  values  Mp  are  required  to  start  up  the  iterative  procedure. 

It  was  observed  that  a  faster  converging  solution  was  achieved  (i.e. 
number  of  iterations  were  reduced),  the  closer  the  guess  values  were  to  the 
actual  solution.  Also  the  user  who  wants  to  use  the  robot,  initially  has 
little  or  no  idea  about  what  guess  values  to  provide.  Therefore  the  user  is 
liable  to  input  guess  values  which  could  lead  to  numerous  iterations  ana 
thereby  waste  expensive  computer  time. 

As  a  remedial  measure,  it  was  decided  that  the  program  instead  or  tne 
user  should  decide  about  the  guess  values.  Thus  the  program  was  modified  into 
a  'Learning  Program'.  The  program  has  a  data  structure,  which  for  different 
initial  to  final  load  positions,  has  the  corresponding  required  end  moment 
values  Mp.  Now  given  an  initial  and  final  load  position,  the  program  searches 
through  its  data  structure  to  locate  the  'guess'  end  moment  values  to  start 
the  iterative  procedure.  Everytime  the  program  is  run,  at  the  end  of  the 
iterations,  the  data  structure  is  updated  appropriately,  thereby  increasing 
the  apparent  knowledge  of  the  program. 

A  FORTRAN  program  has  been  written,  which  carries  out  all  the 
iterative  computations,  matrix  transformations  etc.,  discussed  in  this 
chapter.  The  program  also  takes  care  of  another  important  aspect,  i.e.  the 
user  does  not  have  to  provide  the  initial  guess  end  moment  values  Mp,  to  start 
the  iterative  procedure.  As  discussed  in  the  next  chapter,  the  program 
searches  through  its  database  to  come  up  with  the  requisite  guess  values. 
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The  data  structure  required  to  implement  this  'Learning  Concept', 
could  have  been  appropriately  implemented  in  a  language  like  Pascal  using  tree 
structures.  In  our  case  as  the  remaining  programs  were  in  FORTRAN,  we 
persisted  with  FORTRAN  and  had  to  use  a  five  subscripted  array  as  the  data 
structure. 

Data  Structure  and  Search  Logic: 

As  shown  in  Figure  10,  given  the  coordinates  of  the  initial  and  final 
positions  of  the  load,  we  can  uniquely  define  the  initial  and  final  position 
of  the  load,  by  the  four  parameters  RAD,  $,  0  and  DIS.  The  four  parameters 
are  defined  as  follows 

RAD  =  [xi)2  +  (yi)2]l/2  (20) 

*  =  Tan"1^]  (21) 

yi 

DIS  ■  [ (xf  -  X})2  +  (yf  -  yj)2]1/2  (22) 

l  ( x  ,  -  x . ) 

9  -  COS  h  fms  i~j  (23) 

For  each  set  of  four  such  parameters  we  have  corresponding  end  moment 
values  Mp.  Therefore  when  we  search  for  the  most  accurate  guess  values,  from 
the  available  data  structure,  a  search  through  each  of  the  four  parameters  is 
required.  To  simplify  the  search  procedure,  a  five  subscripted  array 
D( I ,J ,K,L,N) ,  was  used  where  the 

First  subscript  corresponds  to  -  '♦'  values 
Second  subscript  corresponds  to  -  'RAD'  values 
Third  subscript  corresponds  to  -  '0*  values 
Fourth  subscript  corresponds  to  -  'DIS'  values 
Fifth  subscript  corresponds  to  -  end  moment  values  Mp. 
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Now  it  is  shown  in  Figure  11  for  each  '<(>'  we  have  one  or  more 
corresponding  'RAD'  values.  Similarly  for  each  'RAD'  we  have  one  or  more 
corresponding  '0'  values  and  for  each  '0'  we  have  one  or  more  corresponding 
‘ DIS'  values.  Finally  for  each  ‘DIS‘  we  have  the  corresponding  end  moment 
values  for  each  element  of  the  robotic  arm. 

Therefore  given  an  initial  load  position  (x-j ,  y^)  and  a  final  load 
position  (xf,  yf),  using  Equations  (20)  to  (23)  we- can  evaluate  <|>e,  RADe,  0e 
and  DISe  (where  subscript  e  signifies  the  evaluated  parameters  for  which  we 
wish  to  find  best  guess  end  moment  values). 

We  can  now  search  through  the  array  D(I,J,K,L,N).  First,  we  compare 
'  <fre'  with  all  the  '<)>'  values  in  the  array  and  choose  the  one  which  is  closest 
to  '  4>e'  value.  Let  the  chosen  value  be  '<|>c'. 

Next  we  search  through  all  the  'RAD'  values,  corresponding  to  '  4>c *  * 
value,  and  find  the  'RAD'  value  closest  to  'RADe'  value.  Let  this  chosen 
value  be  'RADC'. 

After  this  we  search  through  all  the  '0'  values,  corresponding  to 
‘ RADC '  value,  and  find  the  '0'  value  closest  to  ' 0e '  value.  Let  this  chosen 
value  be  ' 0C‘ . 

Similarly  we  search  through  all  the  'DIS',  corresponding  to  ' 0C ' 
value,  and  find  the  'DIS'  value  closest  to  ' D I Se '  value.  Let  this  chosen 
value  be  'DISC' . 

The  end  moment  values  corresponding  to  ' DISC ' ,  will  be  our  initial 
guess  moment  values  for  the  iterative  procedure.  This  strategy  of  choosing 
the  guess  values  is  seen  to  work  reasonably  well. 

At  the  end  of  the  iterative  procedure,  the  program  computes  the  exact 
end  moments  required  to  move  the  load  from  position  (x^,y^)  to  position 

(xf.yf). 
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The  program  now  updates  its  data  structure  (Array  0( I ,J ,K,L,N) )  and 
4>e»  RADe,  8e  and  DISe  are  inserted  into  the  data  structure  in  the  proper 
places.  As  the  array  has  all  the  values  in  increasing  order  of  magnitude, 
insertion  into  the  array  is  not  difficult.  The  exact  end  moments  computed  by 
the  program  during  the  iterative  procedure  are  also  inserted  into  the  array 
corresponding  to  <j>e,  RADe,  0e  and  DISe. 

After  the  program  has  been  run  for  several  different,  initial  to  final 
load  positions,  the  data  structure  D(I,J,K,L,N) ,  after  being  repeatedly  up¬ 
dated  becomes  sufficiently  diverse.  Now  for  a  varied  range  of  initial  to 
final  load  positions,  the  guess  values  obtained  will  be  sufficiently  accur¬ 
ate. 

The  data  structure  is  implemented  by  having  a  datafile.  Everytime  the 
program  is  run  the  values  for  array  D(I,J,K,L,N)  are  read  in  from  this  data¬ 
file.  At  the  end  of  the  program,  after  the  array  0(I,J,K,L,N)  has  been  up¬ 
dated,  we  overwrite  the  datafile  with  the  present  values  of  array 
U(I,J,K,L,N) .  In  this  way  the  updated  database  is  available  for  the  next  run 
of  the  program. 

If  the  program  has  been  run  repeatedly  which  leads  to  the  data 
structure  becoming  very  large,  the  searching  time  may  be  too  long.  To  avoid 
this,  we  can  easily  edit  the  datafile  and  remove  the  values  in  the  range  which 
are  seldom  used  and  leave  only  those  values  which  are  repeatedly  used.  Also 
after  the  datafile  has  reached  an  optimum  size  and  we  wish  not  to  expand  the 
datafile  anymore,  all  that  Is  required  is  to  remove  the  ‘write  to  datafile 
commands'  from  the  program. 
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NUMERICAL  EXAMPLES  AND  DISCUSSION 

The  program  TYPEL,  listed  in  Appendix  B,  carries  out  the  single 
element  computations  discussed  previously.  For  a  prescribed  end  moment,  the 
end  deflections,  radius  of  curvature  and  position  of  the  neutral  axis  of  the 
deflected  element  are  computed.  This  program  calls  the  IMSL  subroutine  DCADRE 
repeatedly  to  carry  out  numerical  integration. 

The  variable  list  is  as  follows: 


1) 

A1,B2,N  * 

Material  constants  for  nonlinear  material, 

relationship  a  -  Ale  +  B2el/N 

stress  strain 

2) 

81  =  B 

Breadth  of  shim 

3) 

E 

Young's  modulus  of  steel  shim  or  reinforcement 

4) 

H 

Thickness  of  shim 

b) 

HI  =  Hi  = 

Distance  between  neutral  axis  and  point  of 

between  shim  and  tubular  element 

connection 

6) 

HDIS  =  <$h  = 

Horizontal  end  deflection  of  element 

7) 

L 

Length  of  an  element 

8) 

R1  *  Rl  * 

External  radius  of  cylindrical  element 

9) 

R2  =  R2  = 

Internal  radius  of  cylindrical  element 

10) 

ID 

RHO  *  p  = 

VDIS  =  5V  = 

Radius  of  curvature  of  deflected  element 

Vertical  end  deflection  of  the  element 

Vertical  end  deflection  of  the  element 


For  an  element  such  as  shown  in  Figures  1  and  2,  made  of  a  polymeric 
material  with  a  thin  steel  shim,  the  following  numerical  values  were  chosen  as 
input  parameters. 


N  =  0.2 

(dimensionless) 

A1  =  20.0 

psi 

82  =  50,0 

psi 

R1  =  0.525 

i nches 

R2  =  0.475 

inches 

H  =  0.015 

inches 

L  =  2.0 

inches 

B1  «  0.5 

inches 

M  =  3.80 

lb  in. 

E  =  3.0  x  107 

psi 

The  results  of  the  computations  are  as  follows. 


HI  =  0.00749 

inches 

RH0  =  1.549 

inches 

V0IS  =  1.121 

inches 

HDIS  =  0.511 

inches 

It  is  noted  that  the  position  of  the  neutral  axis  (HI)  is  nearly  at  the 
midheight  of  the  steel  shim.  Further,  the  vertical  deflection  at  the  end  of 
the  element  is  about  one-half  the  element  length.  In  this  case  the  internal 
pressure  p  required  to  produce  the  chosen  end  moment  M,  in  terms  of  the 
symbols  of  Figure  2,  is 


M 


3.80 


ttR22(K1  +  Hl)  ir(0.475)2(C525  +  0.00749) 


=  10.1  psi 


Note:  The  remaining  variables  in  the  program  hold  intermediate  values  and  are 
of  no  significance  to  the  user. 


The  program  LEARN,  listed  in  Appendix  B,  carries  out  the  computations 
discussed  in  Chapters  III  and  IV.  It  evaluates  the  end  moments  required  in 
each  of  the  elements  to  shift  the  load  from  one  position  to  another. 

This  program  can  handle  many  different  options.  An  element  in  the  arm 
may  or  may  not  be  pressurized.  The  elements  can  be  of  varying  geometry  and 
the  shim  can  be  at  top  or  at  the  underside  of  the  element.  These  options 
allow  a  great  variation  in  the  motion  of  the  arm. 

Now,  as  previously  discussed,  in  this  program  given  an  applied  end 
moment  on  an  element,  it  was  required  to  know  the  end  deflections  and  radius 
of  curvature  of  the  deflected  element.  To  compute  the  end  deflections  and 
radius  of  curvature  the  previously  discussed  program  TYPEL  is  required.  As 
program  LEARN  goes  through  an  iterative  procedure,  it  would  be  required  to  run 
program  TYPEL  repeatedly.  This  could  be  extremely  time  consuming  and 
expensive.  Instead  we  input  into  the  program  a  two  dimensional  array  A(I,J). 
This  array  consists  of  different  end  moment  values  with  their  corresponding 
end  deflections,  radius  of  curvature,  position  of  neutral  axis,  length  of  the 
element  and  external  radius  of  the  turbular  part  of  the  element. 

The  program  'LEARN'  given  an  end  moment  value  MI,  searches  through 
array  A(I,J)  and  finds  the  two  moment  values  (MA  &  MB)  between  which  MI  lies. 
It  then  linearly  interpolates  between  MA  &  MB  to  evaluate  the  end  deflections, 
radius  of  curvature  and  position  of  neutral  axis  corresponding  to  MI.  This 
explanation  will  become  clear  by  looking  through  the  sample  A(I,J)  input. 


The  variable  list  is  as  follows. 


1)  A(I,J) 


2)  ANG 

3)  D( I ,J ,K,L,M)  = 

4)  H 

5)  HI 

6)  HOIS 

7)  10 1 R 

8)  IMOT 

9)  INCK 

10)  IP1 

11)  IPCK 

12)  IPLT1 

13)  ITYPE 

s 

14)  LN 
lb)  N1 


Array  which  has  for  different  end  moment  values,  the 
corresponding  end  deflections,  radius  of  curvature, 
length  of  element,  position  of  neutral  axis  and  the 
external  radius  of  element 

End  cumulative  rotation  up  to  that  element  in  degrees 
Holds  the  ‘apparent  knowledge'  data  structure 
Thickness  of  shim 

Distance  between  neutral  axis  and  point  of  connection 

between  shim  and  turbular  part 

Horizontal  end  deflection  of  a  particular  element 

Contains  the  IMOT  value  of  previous  element 

-3  means  element  deflects  upwards 

+3  means  element  deflects  downwards 

to  fix  position  of  global  x  axis 

+2  means  element  is  pressurized 

-2  means  element  is  not  pressurized 

0  means  no  pressurized  elements  in  the  arm  up  till  now 

1  means  pressurized  element  has  been  encountered 

-1  means  shim  on  right  side 

+1  means  shim  on  left  side 

1  means  use  array  A(I,J) 

2  means  use  array  B(I,J) 

Length  of  an  element 

Number  of  identical  elements  (geometrical  and  material 
properties  identical  and  also  same  end  moments) 


*, 
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16)  NN  =  number  of  elements  in  the  arm 

17)  NSHIF  =  shift  in  position  of  neutral  axis  between  consecutive 

elements 

18)  NU  =  number  of  elements  with  shim  on  top 

19)  MM  =  number  of  data  sets  in  A(I,J),  1=1  to  MM 

20)  R1  =  outer  radius  of  cylindrical  part 

21)  RHO  =  radius  of  curvature  of  the  neutral  axis  of  an  element 

22)  RL  =  radius  of  load 

23)  SI  =  cumulative  end  rotation 

24)  VDIS  =  vertical  end  deflection  of  a  particular  element 

25)  Ml  =  weight  of  load 

26)  X2,Y2  =  end  point,  values  for  the  current  element  w.r.t.  neutral 

axis 

27)  X3,Y3  =  end  point,  values  for  the  current  element  w.r.t.  mid 

axis  of  element 

28)  XB,YB  =  coordinates  of  c.g.  of  load  at  end  of  each  iteration 

29)  XI, YI  =  coordinates  of  c.g.  of  load  in  its  initial  rest 

position 

30)  XF ,YF  =  coordinates  of  c.g.  of  load  in  the  final  desired 

position 

31)  XORD(I) ,Y0R0(I)=  coordinates  of  end  points  of  the  elements 


Typical  numerical  results  were  obtained  for  line  of  ten  identical 
elements  such  as  shown  in  Figure  4(a)  with  all  the  elements  unpressurized. 

The  requirement  is  to  lift  a  2  lb  load  (wj  =  2)  from  position  (XI,  YI)  = 
(-2.082,  -14.0)  as  shown  in  Figure  4(b),  to  position  (XF,  YF)  =  (8.5,  -8.0)  as 
shown  in  Figure  6(a).  (These  figures  are  not  to  scale.)  The  input  data  is  as 
follows. 


NN 

= 

10 

MM 

= 

10 

NO 

= 

5 

RL 

= 

1.542  inch 

Ml 

= 

2  lbs 

H 

= 

0.015  inch 

XI 

= 

-2.082  inch 

YI 

* 

-14.00  inch 

XF 

= 

8.5  inch 

YF 

s 

-8.0  inch 

The  chosen  dimensions  of  the  five  subscripts  of  array  D(I,J,K,L,M) 

are: 

I  =  1  to  N2  =  3 

J  =  1  to  M2  =  4 

K  =  1  to  N7  =  4 

L  =  1  to  M7  =  4 

M  =  1  to  MX  =13 

The  last  dimension  in  each  subscript  holds  a  dummy  value  (-1000.0),  so  that 
during  the  search  procedure  the  entire  sparse  matrix  need  not  be  searched. 
The  overall  dimensions  were  chosen  to  suit  the  sample  problem  and  not  waste 
memory  space. 
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Specific  data  about  each  element  has  to  be  provided  to  the  program. 

For  example,  IPLT1  specifies  position  of  shim,  IP1  specifies  whether  the 
element  is  to  be  internally  pressurized  or  not.  The  following  table  specifies 
the  properties  of  all  ten  elements  in  this  example  problem. 

I  ■  1  to  NN  NN  =  10  in  this  example 


I  TYPE  (I) 

IP1(I) 

IPLTl(I) 

Nl(  I) 

IMOTl(I) 

1 

+2 

-1 

1 

-3 

1 

+2 

-1 

1 

-3 

1 

+2 

-1 

1 

-3 

1 

+2 

-1 

1 

-3 

1 

+2 

-1 

1 

-3 

1 

+2 

+1 

1 

+3 

1 

+2 

+  1 

1 

+3 

1 

+2 

+1 

1 

+3 

1 

+2 

+1 

1 

+3 

1 

+2 

+  1 

1 

+3 

For  elements  having  similar  geometry  and  stiffness  properties,  we  can 


repeatedly  run  program  TYPEL  to  obtain  the  results  tabulated  below.  Now 
during  the  iterative  procedure  in  program  LEARN,  given  the  end  moment  value, 
we  can  by  interpolation  of  the  array  A(I,J)  obtain  the  corresponding  end 
deflections,  radius  of  curvature  and  the  other  items  in  the  table  below. 


Array  A(I,J) 


I  =  1  to  10 
J  =  1  to  7 


End 

Moment 

(lb-in) 

Radius  of 
Curvature 
p  (in) 

HI 

(in) 

HOIS 

(in) 

Length  of 
VDIS  Element  LN 
(in)  (in) 

Radius  of 
Element  R1 
(in) 

1.02 

5.3979 

0.749255x10-2 

0.0455 

0.3663 

2.0 

0.525 

1.265 

4.3537 

0.749259x10-2 

0.0696 

0.4514 

2.0 

0.525 

1.705 

3.2391 

0.749247x10-2 

0.1247 

0.5981 

2.0 

0.525 

1.915 

2.8904 

0.74924x10-2 

0.1558 

0.6648 

2.0 

0.525 

2.2799 

2.4409 

0.749228x10-2 

0.2164 

0.7745 

2.0 

0.525 

2.4085 

2.3165 

0.749221x10-2 

0.2394 

0.8111 

2.0 

0.525 

2.634 

2.130 

0.7492x10-2 

0.2812 

0.8702 

2.0 

0.525 

3.159 

1.807 

0.74916x10-2 

0.3839 

0.9982 

2.0 

0.525 

3.798 

1.549 

0.7491x10-2 

0.5109 

1.121 

2.0 

0.525 

4.587 

1.346 

0.7489x10-2 

0.6592 

1.232 

2.0 

0.525 

The  numerical  results  are  summarized  as  follows. 


Element  No. 

Total  End  Moment 
to  be  Applied  (lb  in) 

Element 

Pressure  (psi) 

1 

0.181  x 

102 

262.33 

2 

0.173  x 

102 

214.38 

3 

0.149  x 

102 

168.54 

4 

0.114  x 

102 

129.06 

b 

0.721  x 

lQl 

98.86 

6 

0.133  x 

10l 

79.76 

7 

0.502  x 

lQl 

76.24 

8 

0.792  x 

10l 

62.94 

9 

0.792  x 

ioi 

41.96 

10 

0.792  x 

ioi 

20.98 

It  is  to  be  noted 

that  elements  8,  $ 

)  and  10  have  identical  end 

moments.  As  the  exact  load  distribution  is 

not  known. 

these  elements  are 

over-pressurized  so  that 

they  coil  around  and  grip  the 

load  more  tightly  than 

is  required  for  lifting. 

The  pressure  pj  in  the  i-th 

element  whose  end  moment 

is  Mi  was  calculated  from  the  following  recursion  relationships. 

Mi  =  (Pi  -  Pi+l)1^2^  +  Hl) 

It  is  observed  that  the  pressures  in  the  respective  elements  are  higher  for 
the  ones  closer  to  the  support. 

The  elements  in  the  sample  problem  were  identical  but  if  the  arm 
consists  of  elements  such  that  the  cross-sectional  area  of  elements  decreases 
successively  as  they  move  away  from  the  support  (element  1  is  attached  to  the 
support)  then  there  would  be  no  drastic  decrease  in  pressure  between 
successive  elements. 
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There  are  several  possible  extensions  of  this  research.  For  instance, 
the  programming  efficiency  can  be  considerably  improved  if  the  program  is 
rewritten  in  Pascal  which  offers  the  facility  of  tree  type  of  data  structure 
and  dynamic  storage.  By  implementing  a  tree  type  of  datastructure  searching 
and  insertion  of  data  is  greatly  simplified.  Also,  the  mathematical  model  for 
plane  motion  as  discussed  here  can  be  generalized  to  include  out-of-plane 
manipulations.  This  can  be  done  by  the  addition  of  torsion  elements  [6]  at 
the  base  of  element  one.  One  torsion  element  can  move  the  arm  ±180°  out  of 
the  plane  whereas  another  torsion  element  can  move  the  arm  +90°  in  the  plane. 
By  judiciously  pressurizing  these  torsion  elements,  the  robotic  arm  can  have  a 
work  space  consisting  of  volume  of  revolutions  in  each  of  the  four  quadrants. 
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III.  LINEAR  ANALYSIS  OF  UNIFORMLY  STRESSED, 
ORTHOTROPIC  CYLINDRICAL  SHELLS 


ABSTRACT 

Within  the  framework  of  classical  elasticity,  the  nonbuckled  deformations 
are  calculated  for  orthotropic,  right  circular,  thin-walled  cylinders  under 
uniform  load  conditions.  The  principal  direction  of  orthotropy  follows 
parallel  constant  angle  helices.  Nondimensional  system  parameters  involving 
four  material  constants  and  three  loading  conditions  (internal  pressure, 
longitudinal  load  and  pure  torque)  are  identified.  Through  parametric  studies 
deformation  patterns  are  calculated  that  are  unique  to  orthotropy.  Numerical 
examples  illustrate  that  the  proper  selection  of  cylinder  orthotropy  can  lead 
to  designs  with  optimal  deformations  or  load-carrying  capacity.  Results  may 
be  used  for  the  design  of  robotic  actuators  driven  by  internal  pressure. 
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INTRODUCTION 


Improvement  in  the  mechanical  performance  of  cylindrical  or  nearly 
cylindrical  shell -type  structures  may  be  achieved  by  adding  reinforcement  or 
by  making  small  modifications  in  the  basic  structural  geometry.  For  instance, 
the  high  pressure  capacity  of  fire  hoses  is  due  to  the  reinforcing  effect  of 
the  helical  fibers  in  the  rubber  walls.  In  other  cases,  the  buckling 
resistance  to  longitudinal  loading  of  thin-walled  cylindrical  columns  can  be 
greatly  improved  by  the  use  of  a  fluted  design,  or  a  periodic  variation  of  the 
radius  around  the  circumference.  On  the  other  hand,  there  are  cases  where 
high  deformations  are  desirable  as,  for  instance,  in  cylindrical  bellows  used 
as  pipeline  expansion  joints.  Here,  length  changes  of  up  to  100  percent  may 
be  achieved  by  the  use  of  axisymmetric,  periodic  corrugations  along  the  length 
of  a  basic  cylindrical  shape  [1].  Thus,  the  most  flexible  directions  of  the 
bellows  and  the  thin-walled  fluted  column  are  orthogonal. 

The  purpose  of  this  paper  is  to  present  a  unified  continuum  model  of  such 
thin-walled,  cylindrical  shells  using  directional  material  properties.  Such 
continuum  models  are  especially  efficient  for  calculating  overall,  nonbuck¬ 
ling,  elastic  deformations  for  fiber-reinforced  cylinders  [2]  and  for  uniform¬ 
ly  loaded,  anisotropic,  cylindrical  shells  [3,4],  In  the  present  study, 
classical  elasticity  is  employed  to  predict  the  nonbuckled  deformations  for 
uniform,  orthotropic  right  circular  cylinders  with  thin  walls,  subjected  to 
three  types  of  uniform  loads:  a  torque  T  effecting  rotation  about  the  longi¬ 
tudinal  axis,  a  longitudinal  load  P,  and  an  internal  pressure  p.  Strains 
along  the  geometric  axes  of  symmetry  are  presented  as  functions  of  nondimen- 
slonal  system  parameters  involving  these  loads  (applied  separately  and  in  com¬ 
bination),  the  cylinder  geometry,  four  material  properties,  and  the  constant 
helix  angle  80  that  defines  the  direction  of  maximum  stiffness.  See  Figure  1. 


PRINCIPAL  LINE  OF  ORTHOTROPY 


Figure  1 


Definition  of  the  orthotropic  cylinder 


This  continuum  model  may  describe  the  overall  deformation  characteristics  of 
fiber-reinforced  tubes  as  well  as  tubes  with  corrugations  at  arbitrary  helix 
angles,  as  long  as  the  cylindrical  shape  is  maintained  under  load.  Results 
may  be  used  in  future  designs  of  robotic  actuators  that  twist  (40°  <  60  <  80°) 
or  that  only  lengthen  as  for  bellows  (60  *  0)  when  subjected  to  internal 
pressure. 

CONSTITUTIVE  RELATIONSHIPS 

Following  the  development  in  [5],  the  constitutive  law  for  an  ortho¬ 
tropic,  elastic  solid  that  relates  the  six  strain  components  to  the  six  stress 
components  is  given  by 


(1) 


Equation  (1)  is  written  with  respect  to  the  principal  axis  (1,2,3).  The  6x6 
material  matrix,  whose  elements  are  designated  as  a^j,  is  determined  by 
experiments.  The  strains  and  stresses  of  Equation  (1)  are  redefined  using 
single  subscription,  or 


ell  =  el  • 

e22  *  e2 

»  e33  *  e3 

(2a) 

Yi2  85  H  » 

713  *  c5 

»  723  =  e6 

°11  *  °1  • 

®22  =  °2 

,  033  s  °3 

(2b) 

°12  *  °4  • 

°13  *  °5 

,  023  *  °6 

,  Equations 

(1)  become 

ci  “  .  a1  j°j 
J-l 


(3) 


Recognizing  that  the  total  elastic  energy  for  the  orthotropic  solid  is 

invariant  with  respect  to  the  coordinate  system,  the  elements  a'  of  the 

^  J 

material  matrix  in  the  (r,6,z)  coordinate  system  can  be  written  in  terms  of 
a^j  of  the  (1,2,3)  coordinate  system.  That  is 


''J  m*l  n=l 


6  6 

I  I  amn  tyni  ^nj 


(4) 


i  ,j  18  1,2, ...6 

where  q-j j  represents  the  direction  cosines  of  the  (1,2,3)  system  with  the 
(r,8,z)  coordinate  system.  As  shown  in  Figure  1,  the  relationship  is  simply 
one  of  rotation  of  the  (1,2)  axis  to  the  (e,z)  axis  by  an  angle  of  -0O  in  the 
plane  r  «  constant.  The  values  of  q-jj  for  this  rotation  are  listed  in  Table 
1. 


In  the  transformed  coordinate  system,  the  constitutive  law  is 


I 

G 


i 


6 

I 

j*l 


a  o 
ij  J 


which  is  expressed  in  cylindrical  coordinates  as  follows. 


(5) 


71 


For  the  types  of  uniform  loading  considered  here,  the  shear  stresses  o@r  and 
ozr  may  be  approximated  as  zero,  an  assumption  consistent  with  the  elastic 
theory  of  thin-walled  cylinders  [6].  From  Equation  (6),  the  remaining 
stress-strain  relations  are  thus 


The  elements  a  of  the  material  matrix  of  Equation  (7)  are  calculated 
i  j 

from  Equation  (4)  with  a-jj  and  q^-  givea  by  Equation  (1)  and  Table  1, 
respectively.  The  results  are  as  follows*. 


'Note  that  $  in  reference  [5]  is  interpreted  as  -e0  in  Equations 
(8)-(17).  Also,  Equation  (9)  reflects  a  correction  of  the  results  reported  in 
[5]. 
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EVALUATION  OF  MATERIAL  CONSTANTS 


There  are  seven  material  constants  that  appear  in  Equations  (8)-(17), 
namely  E^j,  Eg2»  ^33*  v12*  v23*  V31  ®12*  1^  1 s  necessary  to  determine  how 
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these  constants  are  related  to  each  other  and  then  to  devise  meaningful  tests 
to  measure  them.  To  insure  that  the  strain  energy  is  a  single  valued  function 
of  the  strain  displacements,  the  material  matrix  of  Equation  (1)  must  be 
symmetric.  That  is 


E11  V21  =  E22  v12 

(18) 

e22  v32  =  E33  v23 

(19) 

e33  v13  =  E11  V31 

(20) 

Consider  a  tensile  test  for  which  the  specimen's  longitudinal  axis 
coincides  with  the  longitudinal  or  z-axis  of  an  orthotropic  cylinder,  given 


that  eQ  =  90  deg.  For  a  uniform,  applied  stress  azz  where  a 
it  follows  from  Equations  (7)-(20)  that  the  strains  are 
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Equation  (22)  shows  that  E^  is  the  equivalent  of  Young's  modulus  along  the  1 
axis.  Now  suppose  that  for  60  «  90  deg  the  orthotropic  tube  corresponds  to  a 
fluted,  thin-walled  cylindrical  column  of  a  homogeneous.  Isotropic  material 
with  Young's  modulus  E  and  Poisson's  ratio  v.  This  fluted  column  test  specimen 
thus  exhibits  the  same  stress-strain  behavior  of  the  isotropic  material  in  both 
the  1  (or  z)  direction  and  in  the  3  (or  r)  direction.  The  following 
definitions  of  the  orthotropic  constants  are  consistent  with  Equations 
( 21 ) - ( 23 ) . 


E  "  E11  *  e33 
v  «  v3i  *  v13  «  v12 


(24) 

(25) 


Consider  a  second  tensile  test  for  which  the  specimen's  longitudinal  axis 
again  corresponds  with  the  z-axis  of  an  orthotropic  cylinder,  but  now  let  0O  = 
0  deg.  Again,  the  only  nonzero  stress  is  the  applied,  uniform  tensile  stress 
azz.  It  follows  from  Equations  (7)-(20),  (24)  and  (25)  that 
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Observing  Equation  (27),  it  is  seen  that  the  modulus  of  elasticity  in  the  2  (or 


z)  direction  is  E£2*  Now  let 


E'  =  E, 


v  -  v21  -  v23  (30) 

which  define  respectively  the  effective  Young's  modulus  and  Poisson's  ratio  for 
a  bellows,  consistent  with  the  assumption  that  ogg  *  orr  in  this  tensile  test. 
It  follows,  using  Equations  (24),  (25),  (29),  (30),  that  Equations  ( 18 ) - ( 20 ) 


are  satisfied  if 


v  =  v32  ,  EV  =  E'v 


In  a  third  test,  a  pure  shear  stress  <jqz  is  applied  to  an  orthotropic 
cylinder  for  which  60  *  0.  This  corresponds  to  the  application  of  a  pure 
torque  that  rotates  a  bellows  about  its  longitudinal  axis.  As  for  the  previous 
two  tests,  the  cylinder  is  orthotropic  only  because  of  its  corrugations,  since 
it  is  constructed  of  an  isotropic  material  with  elastic  constants  E,  v  and 
shear  modulus  G.  Torsion  tests  on  such  a  bellows  have  shown  that  [8] 

y6z  s~r  {32) 

which  indicates  that  the  rotation  and  shearing  stress  are  those  predicted  by 
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the  theory  of  thin-walled  tubes  without  corrugations.  When  Equation  (32)  is 
compared  with  the  stress-strain  law  given  by  Equations  (7)-( 17 )  in  this  case, 
it  is  seen  that 

G  =  G12  (33) 

In  summary,  there  are  four  independent  material  constants  needed  to 
describe  the  stress-strain  behavior  of  a  cylinder  with  corrugated  walls  made  of 
an  isotropic  material,  but  modeled  as  an  orthotropic,  smooth-walled  tube.  The 
stiffest  direction  (modulus  E)  follows  the  corrugation  lines  forming  constant 
angle  helices.  Orthogonal  to  the  helices  is  the  weakest  direction  (modulus 
E').  The  independent  material  constants  are  E,  E',  G  and  v.  Given  uniform 
stresses  o6e,  o2Z,  orr,  and  o9z,  the  uniform  strains  are  calculated  from 
Equation  (7).  Using  Equations  (24),  (25),  (29-31)  and  (33),  Equations  (8)-(17) 
give  the  elements  of  the  material  matrix  as  follows. 
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These  results  can  be  used  to  calculate  the  strains,  displacements  and  rotations 
of  the  cylinder  for  the  special  cases  of  loading  discussed  below. 


EQUILIBRIUM  AND  COMPATIBILITY 


As  shown  in  Figure  1,  the  thin-walled  cylinder  is  subjected  to  three  types 
of  uniform  loads:  an  internal  pressure  p,  a  longitudinal  load  P  centered  on 
the  z-axis,  and  a  pure  torque  T  causing  rotation  about  the  z-axis.  Since  there 
are  no  boundary  or  edge  constraints,  the  resulting  uniform  stresses  may  be 
derived  using  equilibrium  conditions  and  elementary  methods.  In  terms  of  the 
mean  radius  R0  and  shell  thickness  t0,  which  are  essentially  constant  during 
loading,  the  stresses  are: 
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Consistent  with  the  thin-wall  assumption,  the  radial  stress  is  negligible,  or 
orr  s  0  [6]. 

Since  these  stresses  and  their  corresponding  strains,  as  given  by 
Equations  (7)  and  (34-43),  are  all  uniform,  they  are  all  independent  of  the 
cylindrical  coordinates  (r,6,z).  Thus  the  81  St.  Venant  strain  compatibility 
conditions  are  automatically  satisfied  [7], 
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STRAIN-DISPLACEMENT  RELATIONS 


The  general  strain-displacement  equations  in  polar-cylindrical  coordinates 
given  in  [7]  are  simplified  as  follows  to  express  the  condition  that  there  are 
no  variations  in  displacement  in  the  0  or  circumferential  direction. 
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Here,  Ur,  U0,  and  U2  are  respectively  the  displacements  in  the  r,  e  and  z 


di rections. 


Following  a  procedure  in  [3],  these  relationships  can  be  integrated  in 
terms  of  three  functions  fi(r,z),  f£( z )  and  f3(z)  which  are  to  be  determined. 
From  Equation  (47): 

Ur  =  r  err  +  fj(r,z)  (51) 

It  follows  from  Equations  (48)  and  (51)  that: 

f i ( r ,z )  *  r(e00  -  epr)  (52) 

From  the  last  two  equations,  then, 

Ur  *  r  Eee  =  ^oc60  (53) 

where  r  *  R0  for  a  thin-walled  cylinder. 

From  Equation  (49): 

Uz  *  z  eZ2  +  f2(z)  (54) 

where  the  second  function  of  integration,  f2(z),  is  Independent  of  r  since  the 
cylinder  wall  is  very  thin  and  any  variations  across  this  thickness  have  no 
significant  effect  on  the  longitudinal  deformation.  Further,  by  imposing  the 
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condition  of  no  rigid  body  motion,  where  Uz  =  0  at  z  *  0,  and  noting  that  Uz 
can  be  at  most  a  linear  function  of  z,  then  it  is  apparent  that  f2 ( z )  *  0. 

Thus  Equation  (54)  thus  yields  the  position  5  =  z  +  Uz  of  a  material  point 
originally  at  point  z  measured  from  one  end  of  the  cylinder,  or 

i  -  (1  +  tzz)z  (55) 

The  circumferential  displacement  is  obtained  by  integrating  Equation  (50) 
and  then  using  Equation  (55).  That  is 

Ue  '  Hez  +  f3(z)  (56) 

Here  the  function  of  integration  is  again  independent  of  r,  which  is  consistent 
with  the  thin-wall  assumption.  By  imposing  the  condition  of  no  rigid  body 
rotation  where  U@  *  0  at  z  *  0,  and  noting  that  U0  can  at  most  be  a  linear 
function  in  z,  it  is  evident  that  f3(z)  *  0.  Thus  Ue  *  6  y0z. 

In  summary,  the  radial,  circumferential  and  longitudinal  displacements  are 
given  in  terms  of  the  strains,  the  mean  radius,  and  the  z-coordinate  only,  or 
ur  *  Ro  e08  *»  U0  *  (1  +  ezz)zY0z  ;  Uz  =  z  tzz  (57) 

It  is  noted  that  Uq  and  Uz  are  the  same  displacements  assumed  [2]  for  the  solid 
orthotropic  cylinder  rotating  about  the  z-axis. 

A  further  displacement  of  practical  interest  is  the  angle  of  rotation  $ 
for  a  cross  section  of  the  cylinder  at  the  coordinate  z.  Within  the  limits  of 
linear  theory,  $  =  U0/Ro.  Thus,  from  Equation  (57): 

♦  ■  f  d  *  *22>*ez  (58) 

0 

PARAMETRIC  STUDIES 

All  of  the  calculations  for  the  parametric  studies  that  follow  are  based 
on  Equation  (7),  the  stresses  of  Equations  (44)-(46),  and  the  elements  of  the 
material  matrix  given  by  Equations  (34)-(43).  The  specific  values  of  E/G  «  3 
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and  v  *  0.5  were  chosen  since  they  are  characteristic  of  the  polymeric  and 
rubber-like  materials  to  be  used  in  future  applications. 


Effects  of  Loading  without  Constraint.  In  the  first  series  of  parametric 
studies,  the  separate  effects  of  each  loading,  p,  P  and  T  on  the  homogeneous 
strains  Yez  and  e22  were  investigated  for  an  unconstrained  cylinder.  The  goal 
was  to  determine  peak  strains  as  a  function  of  the  appropriate  nondimensional 
system  parameters. 

The  three  nondimensional  parameters  involving  the  shear  strain  Yez  are 
defined  by  the  ordinates  of  Figures  2  through  4,  in  which  the  only  nonzero  load 
is  p,  P  and  T,  respectively.  For  the  isotropic  case  (E/E '  =  1),  all  three  of 
these  figures  show  no  shear  strain  variations  with  the  helix  angle  60,  a  result 
that  could  be  anticipated.  For  p  loading  only  (Figure  2)  and  for  P  loading 
only  (Figure  3)  the  negative  ordinates  indicate  that  Yez  is  negative,  or  that 
the  cylinder  "unwinds"  because  of  orthotropy  (E/E 1  >  1).  For  T  loading  only 
(Figure  4),  the  positive  ordinate  indicates  that  the  cylinder  "winds  up"  if  T 
is  applied  as  shown  in  Figure  1.  The  interesting  result  is  that  the  peak  value 
of  each  respective  shear  strain  parameter  occurs  at  a  distinctly  different 
value  of  e0.  That  is: 

E  t  Yg 

peak  ( — p — — )  occurs  for  6Q  *  53  deg 


E  t  R  Yg 

peak  ( - p-2 — — )  occurs  for  6Q  «=  30  deg 


E  t  R  Yg 

peak  ( - 2-^ - -)  occurs  for  0Q  *  45  deg 

For  each  load  type,  there  is  no  shift  in  these  respective  values  of  90  for 
1  <  E/E'  <  100. 
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HELIX  ANGLE  B0,  DEG. 


Figure  3  Unrestrained  twisting  due  to  longitudinal  load  (left  ordinate) 

and  unrestrained  longitudinal  strain  due  to  pure  torque  (right 
ordinate) 
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The  three  nondimensional  parameters  involving  the  longitudinal  strain  e22 
are  defined  by  the  ordinates  of  Figures  5,  6  and  3  (right  ordinate),  given  that 
p,  P  and  T  are  the  only  nonzero  loads,  respectively.  For  the  isotropic  case 
( E/E  *  c  1),  e22  is  zero  for  two  load  cases:  with  the  p  load  as  a  direct 
consequence  of  choosing  v  as  0.5,  and  with  the  T  load  as  expected  from 
elementary  solutions.  With  orthotropy,  the  peak  values  of  the  longitudinal 
strain  parameters  are  as  follows: 

E  t  e 

peak  ( — "°q— • ■■■)  occurs  for  e0  =  0 
^  o 


E  t  R  e 

peak  ( - — — )  occurs  for  S£ 


peak  (-  R°  *2Z 


)  occurs  for  6  *  30  deg 


T  ' -  —  o 

From  Figures  5  and  6  it  is  observed  that  the  values  of  these  parameters  are 
only  somewhat  depressed  for  0  <  80  <  15  deg,  but  that  this  depression  increases 
more  rapidly  as  E/E 1  increases.  From  Figure  3  (right  ordinate)  the 
longitudinal  strain  is  seen  to  be  negative,  indicating  the  tendency  of  the 
orthotropic  cylinder  to  shorten  as  it  winds  up  under  a  positive  torsion  load. 

Effects  of  Radial  and  End  Constraint.  While  the  studies  above  dealt  with 
strain  behavior  for  loadings  applied  one  at  a  time,  one  can  imagine  a  multitude 
of  practical  cases  involving  combined  loads.  Consider  now  one  such  case  where 
an  orthotropic  cylinder  is  required  to  operate  in  close  proximity  with  other 
mechanical  parts.  To  avoid  longitudinal  and  radial  or  circumferential 
expansions  when  pressurized,  suppose  that  the  unloaded  orthotropic  cylinder 
just  fits  within  the  confines  of  a  rigid,  closed-end  tube.  Under  loading, 
then: 


ezz  "  e96  B  0 


(59) 
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The  confining  tube  has  frictionless  walls  where  its  design  allows  for  the 
application  of  a  torque  T  to  the  orthotropic  cylinder.  The  confining  tube  does 
not  inhibit  torsion.  From  Equations  (7)  and  (44-46),  the  two  constraints  of 
Equation  (59)  become: 


.1*^0  •  P  T 

+  2  a22  t  +  a22  2*R  t  a26  - D2+ 

0  0  2irRot0 


/  '  .1  1  \  ^  .  < 
(all  +  2  a12)  t  +  a12 
o 


2*V  +  aie  zJt 
0  0 


The  compatible  equation  relating  the  loadings  p  and  T  is  found  by  eliminating  P 
between  Equations  (60)  and  (61),  or 


2,  (a12  *12  ~.ail  ^2,  .  s 
a16  a22  *  a12  a26 


(62) 


The  pressure-longitudinal  load  parameter  is  found  by  eliminating  T  between 
Equations  (60)  and  (61).  In  terms  of  S  of  Equation  (62),  this  is 


,a12  .  1  *26 


i?  -  -  (-- +  i +  s) 


The  results  of  Equations  (62)  and  (63)  are  shown  in  Figures  7  and  8. 


Figure  7  shows  that  the  pressure-torque  parameter,  the  inverse  of  Equation 
(62),  has  a  peak  value  at  80  =  50  deg  for  all  E/E '  >  1.  For  an  imposed  torque, 
this  gives  the  maximum  required  cylinder  pressure  for  full  cylinder 
confinement.  Figure  8  (the  right  ordinate)  shows  that  the  end  load-pressure 
parameter  defined  by  Equation  (63)  is  independent  of  E/E'  if  that  ratio  exceeds 
unity.  For  a  fixed  pressure  level  p,  the  compatible  end  load  P  always 
decreases  with  0O  to  maintain  full  cylinder  confinement.  It  is  noted  that  the 
isotropic  cylinder  is  not  shown  in  Figure  8  (the  right  ordinate)  since  the  only 
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PRESSURE -ROTATION  PARAMETER, 


END  LOAD -PRESSURE  PARAMETER, 


way  that  full  cylinder  confinement  can  be  met  for  E/E  *  1  is  that  P  «  p  *  0. 
For  this  latter  set  of  conditions,  it  is  recalled  that  a  nonzero  value  for  T 
always  implies  full  cylinder  confinement;  that  is.  Equation  (59)  is  satisfied. 

Using  the  results  of  Equations  (62)  and  (63)  it  is  possible  to  calculate 
the  helix  angle  that  will  minimize  the  rotation  +  (minimize  Yez)  at  a  given 
pressure  level.  The  conclusion  in  this  case,  based  on  the  results  of  Figure  8, 
left  ordinate,  is  that  60  “  53  deg  will  produce  a  minimum  $  for  all  E/E'  >  1. 

Effects  of  Torsional  and  End  Constraint.  Consider  now  the  special  case 
where  the  orthotropic  cylinder  is  fully  restrained  from  rotation  and  also  from 
longitudinal  displacement.  That  is 

Y6z  *  Ezz  "  0  (64) 

With  Equations  (7)  and  ( 34 ) - ( 46 )  these  two  constraint  equations  are: 


,  .  l*,pRo  •  P  •  T 
(a12  +  2  a22}  t  +  a22  2wR  t  +  a26  „_n2fc  *  0 


2*R  t 
oo 


•  .  1  <  ,  p  o  •  P  •  T 

l16  2  a26}  t  +  a26  2wR  t  +  a66  TTT” 

o  00 

o  o 


By  eliminating  P  between  Equations  (65)  and  (66),  the  torque-pressure  parameter 


is  deduced  as: 


-2.  •?6  '  I6  ‘I2)  •  S’ 

a22  a66  ‘  a26  a26 


The  compatible  end  load-pressure  parameter  is  found  by  eliminating  T  between 
the  same  two  equations.  In  terms  of  S'  .of  Equation  (67),  this  is 


-2.  (Iliads') 


The  results  of  these  studies  are  shown  In  Figures  9  and  10.  For  a 
constant  pressure  p  and  a  fixed  modulus  ratio  E/E',  it  is  observed  that  both 
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the  applied  torque  T  and  end  load  P  reach  peaks  that  are  highly  dependent  on 
the  helix  angle.  For  Instance,  If  the  orthotropic  cylinder  were  to  be  used  as 
an  efficient  torque -exerting  device,  we  would  pick  60  *  60  deg  if  E/E  *  *  2,  but 
would  pick  0O  «  80  deg  if  E/E  *  *  50.  However,  If  this  cylinder  were  to  be  used 
as  a  longitudinal  or  end  loading  device,  we  would  choose  0O  *  62  deg  for  E/E '  - 
50.  These  results  are  contingent,  of  course,  on  the  requirement  that  all 
rotations  and  longitudinal  motions  be  completely  supressed. 

There  are  three  observations  for  this  case  where  0O  *  90  deg  and  p  >  0. 
First,  no  end  load  is  required  for  the  longitudinal  strain  to  vanish,  which  is 
a  direct  consequence  of  having  chosen  v  »  0.5.  Second,  no  torque  is  required 
and  fez  automatically  zero,  a  result  that  agrees  with  intuition.  Third,  a 
straightforward  calculation  of  the  circumferential  strain  reveals  that  e0e 

reaches  Its  peak  value  when  0O  *  90  deg,  for  all  E/E'  >  1. 

# 

DESIGN  EXAMPLES 

The  following  examples  show  the  utility  of  the  above  parametric  studies  in 
the  design  of  special  purpose  configurations.  The  finite  length  cylinders  are 
assumed  to  be  fitted  with  end  caps  that  have  negligible  radial  stiffness  so 
that  homogeneous  stresses  are  maintained  throughout.  The  material  and 
geometric  properties  common  to  all  three  design  examples  discussed  below  are  as 
follows. 

E  *  1.5  x  107  N/m2  ,  E'  -  1.5  x  106  N/m2, 

G  ■  5  x  10®  N/m2  ,  v  ■  0.5, 

t0  *  1.5  mm  ,  R0  ■  15  um  ,  lQ  •  50  mn. 

In  the  first  design,  0O  ■  30  deg  and  the  maximum  allowable  Internal 
pressure  Is  p  ■  1.9  x  10^  N/m2.  What  Is  the  maximum  torque  T  that  such  a 
cylinder  can  exert  to  loosen  a  bolt,  for  Instance?  With  E/E'  ■  10,  the  value 
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of  the  shear  strain-pressure  parameter  for  unconstrained  rotation  (T  «  0)  is 
given  in  Figure  2  as  -4.87,  from  which  Yez  *  -0.62  rad.  The  maximum  torque 
capacity  is  then  the  value  of  T  for  which  yqz  *  +0.62  rad.,  corresponding  to 
complete  supression  of  rotation.  With  this  shear  strain  and  the  ordinate  of 
--E4 gu re-4~t— for  wh  1  Ch~th~e  shear  strain-torque  parameter  is  1.55  at  0O  *  30  deg 
and  E/E1  *  10,  the  maximum  possible  torque  is  thus  calculated  as  T  *  2.02  N-m. 

For  this  first  design,  what  is  the  length  change  in  the  cylinder  for  this 
maximum  torque  condition?  From  Figure  5,  the  longitudinal  strain-pressure 
parameter  is  4.2,  for  which  ezz  ■  0.53  for  unconstrained  rotation  (T  *  0).  For 
fully  constrained  rotation  (Yez  *  0)  T  *  2.02  N-m.  The  corresponding  value 
of  longitudinal  strain  from  the  right  ordinate  of  Figure  3  is  calculated,  or 
czz  *  -0.37.  By  superposition,  the  net  result  Is  ezz  *  0.16,  the  length  change 
of  the  cylinder  is  Uz  *  8  nm  and  the  final  cylinder  length  is  6  *  58  mm,  which 
are  found  from  Equations  (57)  and  (55)  respectively  for  z  *  t0  *  50  nm. 

In  the  second  design,  the  maximum  allowable  Internal  pressure  is  again 
p  *  1.9  x  105  N/m2.  What  value  of  60  will  produce  the  maximum  end  load  P?  The 
procedure  is  to  allow  the  cylinder  to  expand  freely  in  the  longitudinal 
direction  with  a  strain  eZ2  compatible  with  Figure  5;  and  then  for  the  same  0O 
to  apply  a  compressive  load  P  compatible  with  Figure  6  where  this  load 
completely  nullifies  the  previous  strain  ezz.  This  is  equivalent  to  finding 
the  value  of  0O  that  maximizes  the  ratio  Y5/Y5,  or  the  ratio  of  ordinates  of 
these  two  respective  Figures.  In  these  terms,  the  peak  compressive  value  of  P 


is  given  by 


Yc  o 

P  -  Oy2*)  PR«  «  341  N 
6  max 


corresponding  to  80  «  58  deg. 
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In  the  third  design,  the  applied  torque  and  the  longitudinal  strain  are 
zero.  What  value  of  60  will  produce  the  maximum  cylinder  rotation  at  an 
arbitrary  internal  pressure?  In  this  case,  Yez  the  sum  oY  two  components: 
that  due  to  the  internal  pressure  (Figure  2)  and  that  due  to  the  longitudinal 
load  (Figure  3,  left  ordinate).  In  terms  of  the  ordinates  Y£  and  Y3  of  these 
respective  Figures,  as  well  as  the  ordinates  Y5  and  Yg,  the  total  shear  strain 
reduces  to: 

Y,  P 

Tez  *  <-Y2  +  rr  (70) 

60 

In  Equation  (70),  the  value  of  P  was  eliminated  from  Y3  by  the  multiple  Yg/Yg. 
Using  a  trial  and  error  procedure,  the  peak  value  of  the  bracket  term  in 
Equation  (70)  was  calculated  as  -4.94,  for  which  60  *  67  deg.  If  the  internal 
pressure  is  p  «  3.8  x  104  N/m*,  then  the  corresponding  shear  strain  is  Yez  * 
-0.125  rad;  and  the  rotational  displacement  and  angle  of  rotation  at  z  =  t0, 
calculated  from  Equations  (57)  and  (58),  are  Ug  *  -6.25  irm  and  ♦  «  -0.42  rad, 
respectively. 


SUMMARY  AND  CONCLUSIONS 


Within  the  limits  of  linear  theory,  the  deformation  behavior  of  uniformly 
stressed,  orthotropic  cylindrical  shells  is  described  by  several  nondimensional 
parameters  involving  four  Independent  material  constants  E,  E',  G,  v,  and  three 
loads:  Internal  pressure,  longitudinal  load  and  pure  torque.  It  is  assumed 
that  buckling  is  absent  and  that  the  cylinder  is  sufficiently  constrained  so 
that  the  cylindrical  shape  is  always  maintained.  As  shown  in  Figure  1,  the 
principal  directions  of  orthotropy  follow  constant  angle  helices.  The  more 
Important  findings  of  the  parametric  studies,  unique  to  this  type  of 
orthotropy,  are  summarized  as  follows. 


1.  Longitudinal  loading  P  and/or  internal  pressure  p  effects  twisting 
about  the  longitudinal  axis. 

2.  For  P  loading  only,  maximum  twisting  deformation  occurs  for  an 
orthotropy  angle  of  60  ■  30  deg;  but  for  p  loading  only,  this  occurs 
when  e0  *  53  deg. 

3.  For  a  pure  torque  T,  maximum  twisting  deformation  occurs  when  0O  =  45 
deg. 

4.  A  pure  torque  T  produces  a  change  in  cylinder  length,  which  is  not  the 
case  for  isotropic  cylinders. 

5.  The  cylindrical  strains  Yg2,  699  and  e22  are  all  amplified  for  a  given 
set  of  loads  (p,  P,  T)  as  the  magnitude  of  the  orthotropy  increases, 
that  is,  as  E/E '  increases,  given  that  the  geometry  and  other  material 


constants  remain  fixed. 

6.  The  maximum  longitudinal  strain  e22  for  P  and  p  loading  occurs  for  eQ 
■  0;  but  e2Z  is  maximum  for  T  loading  when  60  =  30  deg. 

7.  The  design  examples  show  that  the  proper  selection  of  cylinder 
orthotropy  (80  and  E/E ' )  can  lead  to  optimal  deformations  or 


load-carrying  capacity.  Such  designs,  however,  are  highly  sensitive 
to  both  the  type  of  loading  and  strain  constraint  conditions  . 

The  results  of  these  studies  are  important  to  the  design  of  versatile  and 
efficient  pressure-controlled  actuators  made  of  rubber-like  materials  that  can 
sustain  relatively  high  strains.  Orthotropy  may  be  achieved  by  corrugations  of 
the  wall.  Robotics  is  a  logical  field  of.  application  for  such  actuators. 
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NOMENCLATURE 


E 

E* 

Ell 

6 

6ij 

P 

Ro 

T 

Yk 

aij 

fl(r,z) 
f2(z).  f3(z) 
Ao 
P 

r,0»z 

*o 

Ur 

Uz 

Ue 

5 


Young's  modulus  in  strong  direction  of  orthotropy 

Young's  modulus  in  weak  direction  of  orthotropy 

Modulus  of  elasticity  in  i  direction  where  i  is  one  of  the 
orthotropic  axes,  i  »  1,2,3 

Shear  modulus  of  an  isotropic  material  with  Young's  modulus  E 

Shear  modulus  in  1-j  plane,  where  i  and  j  are  orthotropic  axes, 
(i,j)  -  1,2,3,  i  *  j 

Longitudinal  load 

Mean  radius  of  cylinder 

Applied  end  torque 

Value  of  ordinate  for  Figure  k,  k  «  2,3, ...6 
elastic  constants  In  orthotropic  coordinate  system 

(1.J)  »  1,2, ...6 

Elastic  constants  In  cylinder  coordinate  system 

(1 , j )  *  1,2, • . .6 

Function  of  Integration 

Functions  of  integration 

Length  of  cylinder 

Internal  pressure 

Direction  cosine  terms  relating  the  (1,2,3)  axis  to  (r,e,z) 
axis,  (i,j)  *  1,2, ...6 

Cylindrical  coordinates 

Wall  thickness 

Radial  displacement 

Longitudinal  displacement 

Circumferential  displacement 

Longitudinal  coordinate  of  material  point  where  5  a  z  at 


Shear  strains  In  i-j  plane  in  orthotropic  coordinate  system, 
(i,1)  -  1,2,3,  i  *  j 

Shear  strain  in  z-r  plane 

Shear  strain  in  8-r  plane 

Shear  strain  in  8-z  plane 

Strains  in  orthotropic  coordinate  system,  single  index 
notation  i*l,...6 

Strains  in  cylinder  coordinate  system,  single  index  notation 
i*l,...6 

Strains  in  orthotropic  coordinate  system,  double  index 
notation  1  *  1,2,3 

Strain  in  radial  direction 

Strain  in  longitudinal  direction 

Strain  in  circumferential  direction 

Angle  between  strong  axis  of  orthotropy  and  circumferential 
direction  on  cylinder 

Poisson's  ratio  of  isotropic  material  with  Young's  modulus  E 

Ratio  of  extensional  strain  in  i  direction  to  extensional 
strain  in  j  direction  for  orthotropic  coordinate  system, 
(1.J)  •  1.2,3,  i  *  j 

Stresses  in  orthotropic  coordinate  system,  i*l,2,...6 
Stresses  in  cylinder  coordinate  system,  i  «  1,2,... 6 

Stresses  in  orthotropic  coordinate  system,  (i,j)  *  1,2,3 

Normal  stress  in  radial  direction 

Shear  stress  in  z-r  plane 

Normal  stress  in  longitudinal  direction 

Shear  stress  in  8-r  plane 

Shear  stress  in  s-z  plane 

Normal  stress  in  circumferential  direction 

Angle  of  rotation  of  cross-section 


REFERENCES 


1.  Wilson,  J.  F.,  “Mechanics  of  Bellows:  A  Critical  Review,"  International 
Journal  of  Mechanical  Sciences,  1985  (in  press). 

2.  Verma,  P.  0.  S.  and  Rana,  0.  H.,  "Rotation  of  a  Circular  Cylindrical  Tube 
Reinforced  by  Fibres  Lying  Along  Helices,"  Mechanics  of  Materials,  Vol.  2, 
pp.  353-359,  1983. 

3.  Reissner,  E.,  "On  Uniform  Stress  and  Strain  in  Axially  Homogeneous 
Cylindrical  Shells,"  International  Journal  of  Solids  and  Structures,  Vol. 

6,  pp.  133-138,  1970. 

4.  Reissner,  E.  and  Tsai,  W.  T.,  "Pure  Bending,  Stretching,  and  Twisting  of 
Anisotropic  Cylindrical  Shells,"  Journal  of  Applied  Mechanics,  Vol.  41,  pp. 
168-174,  1974. 

5.  Lekhnitskii,  S.  G.t  Theory  of  Elasticity  of  an  Anisotropic  Elastic  Body, 
Translated  from  Russian  by  FT  Fern,  Holden-Day,  Inc.,  San  Francisco,  1963. 

5.  Timoshenko,  S.  P.  and  Woinowsky-Krieger,  S.,  Theory  of  Plates  and  Shells, 
2nd  Edition,  McGraw-Hill,  New  York,  1959. 

7.  Sokolnikoff,  I.  S.,  Mathematical  Theory  of  Elasticity,  2nd  Edition, 
McGraw-Hill,  New  YorF,"  T9K. - - 

.  Dahl,  N.  C.,  "Toroidal -Shell  Expansion  Joints,"  Journal  of  Applied 
Mechanics,  Vol.  20,  pp.  497-503,  1953.  - 


( 


100 


* 


APPENDIX  A 

COMPUTATION  OF  NEUTRAL  AXIS 


The  iteration  steps  to  compute  the  position  of  the  neutral  axis  of  the 
composite  section  of  an  element,  given  an  end  moment  Ma,  were  briefly  outlined 
in  Chapter  2.  The  iteration  steps  including  various  formulae  derivations,  are 
now  discussed  in  detail. 

To  set  up  the  force  equilibrium  equation  for  the  composite  section  it 
is  needed  to  integrate  over  elemental  areas  across  the  section.  For  the 
cylindrical  part  as  shown  in  Figure  2a,  the  area  of  elemental  steps  are  as 


follows. 


where 


area  of  strip  ABCD  is  2 J  Ari  dy 

area  of  strip  EFGH  is  J  Ari  dy  - J  Ar2  dy, 

Arl  ■  Rl2  -  (Rl  -  y  +  Hi)2 
Ar2  =  R22  -  (Rl  -  y  +  Hl)2 


(A-l) 

(A*2) 


The  iteration  steps  are 


axis. 


a.  Choose  an  initial  value  of  Hi  to  fix  the  position  of  the  neutral 


b.  After  the  neutral  axis  position  is  chosen  the  radius  of  curvature 


p  of  the  deflected  element  is  computed.  Due  to  applied  end  moment,  the  part 
of  the  section  above  the  neutral  axis  will  be  in  tension,  whereas  below  the 
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neutral  axis  the  section  will  be  in  compression.  Refer  to  Figures  2a  dnd  2c. 


For  the  linear  reinforcement  material  the  stress  strain  relationship  is 


a  =  Ee 


and  strain  curvature  relationship  is 


e  =  y/p 

From  Equations  (A. 3)  and  (A. 4)  we  get 

o  =  E  y/p 

For  the  nonlinear  synthetic  material  the  stress  strain  relationship  is 

o  *  Ae  +  Bel/n 

From  Equations  (A.4)  and  (A.6)  we  get 

•  -  A#  ♦  B(^)1/n 


(A. 3) 


(A.4) 


(A. 5) 


(A.6) 


(A. 7) 


F  *  Force  *  /  odA 
A 


(A. 8) 


From  Equations  (A. 5)  and  (A.8),  the  compressive  force  in  the  lower  part  of  the 
reinforcement  (below  the  neutral  axis)  is 


o  -Eb(H  -  HjT 

F  =  /  E  b  dy  = - 5 - - - 

H-Hj  p  2p 

Similarly  the  tensile  force  in  the  upper  part  of  the  reinforcement  is 


(A. 9) 


Hfl  y  EB  H‘ 

F  •  J1  E  *  b  dy  ■  -jji 


(A. 10) 


From  Equations  (A. 7)  and  (A.8),  the  tensile  force  in  the  cylindrical  part  is 


f  =  —  (f)  +  B(i)1/n(I) 


(A. 11) 


Y  =  2AC  J 

H, 


^1”**2+R1 


*»*  l  „  iU  -F7i^ 


Rj-Rg+Hj 


2R1+H1 


r1+r2+h1 


*7V I  dy] 


(A. 12) 


R, -R?+H 

i  -  a;  1  y1/n 

Hi 


R1+VH1 


»*'  TTi  * +  i  „  iU »  {y*rt  -  a?* 


Rj-Rg+Hj 


2R1+H1 


Rl+R2+Hl 


y1/n/^T  <w 


(A. 13) 


Y  and  1  are  evaluated  using  numerical  integration  techniques. 


To  have  force  equilibrium  we  know  that  the  sum  of  compressive  and 


tensile  forces  is  equal  to  zero.  Thus  from  Equations  (A.9),  (A. 10)  and  (A.ll) 


we  get 


Eb(H  -  H,)2  Eb  H  2  ,  ,  1/n 

■— irJ“  +  -rr-  +  iY  +  B^  / 


Eb(H4  -  2HH.)  -  2y  , 

- inr - B(i)1/n  I 


2  I.  i 

Eb(H  -2HH. )  -  2y  in1 

- mr - <?) 


1  _  rEb(H2-2HH,)  -  2t,  £ 
p  L  2b!  j 


r  2B1  i 

P  =  L - 2 - j 

Eb(H^-2HH1)-2Y 


(A. 14) 


>-VV-Y-,“-V-  ‘ - V-  > '» V- Y-Y • V V 


c.  For  the  computed  radius  of  curvature  p  and  chosen  position  of 

neutral  axis  the  moment  of  resistance  of  the  section  Mc  is  evaluated. 

Moment  of  resistance  *  /  ay  dA  (A. 15) 

A 

Using  Equations  (A. 5)  and  (A. 15)  the  moment  resisted  by  the  reinforcement  is 

M  -  I1  y2  bdy  +  /  y2  I  bdy 
o  p  H1~ H  p 

Mp  =■£(<*)  (A. 16) 

where 

“  *  f^-Hl3  +  (H  "  Hj)3]  (A. 17) 

Using  Equations  (A. 7)  and  (A. 15)  the  moment  resisted  by  the  nonlinear 
cylindrical  part  is 

*  <i>  *  +  <i>l/n''  (A. 18) 

where 


and 


R.-R-+H.  , _  n 

6  *  2AC/1  2  1  y2/TT  dy  +  f1 

H  «  rl  o 


R.+Rg+Hj 


1 


2R.+H. 

+  /  1  1 

Rl+R2+Hl 


VVhi 


y  f\i dy] 


(A. 19) 


Ri-Ro+Hi  —  R.+R^+H, 

n  *  2BC/1  2  1  y  n  /XT  dy  +  f1  2  1 

v  1  R1-R2+H1 


n+1 

n 


H, 


2R.+H, 

*  j1  1 

Rl+R2+Hl 


n+1 


(A. 20) 


The  total  moment  resisted  by  the  section 


up 


Mc  -  Mr  +  Mn 

From  Equations  (A. 16)  and  (A. 18)  we  obtain 

Mc  =  (i)  a  +  (i)  8  +  (i)1/n  n  (A. 21) 

d.  If  the  following  inequality  holds  for  a  specified  tolerance  limit 
e 

M  -  M 

_a _ £  < 

"a 

then  the  iterations  stop.  Otherwise  a  new  position  of  neutral  axis  is  chosen 
and  steps  (b)  to  (d)  are  repeated. 


s 

fm 

-> '->V 


APPENDIX  B:  PROGRAM  LISTINGS 


PROGRAM  TYPE! 


//  JOB 

//  EXEC  FTVCG 

T/c.sysxn  66  • . 

c 


REAL'S  HI. H.N.B1.M.A1.B2.R1.R2. RHO. RHIN, ALPHA. 
1BETA, NETA.KSITA. GAMA. Ml. M2. U.Z. HI. L. MM. 1NTGRL. 
2E, VOIS. HOIS. CHK.O.R 
re a6(  i .  io)  n.b i . h . 6 i  ,1a '."t .  a i . B2 .m . i . 

10  FORMAT (  F  3 . 1 , F 4 . 1 , F7 . 4 , 2F6 . 3. F4  . 1 . F5 . 1 . F5 . 1 . ES . 2. 
_  WRXTE__( 3,10)  N.B1.H.R1.R2.L.A1.B2.M.E 

"c  . . . .  . . ; .  . . . 

C  LOOP  TO  VARY  THE  APPLIED  MOMENT 

C  __ 

“do  ao  '1-1. T  .  . . . . .  . 

WRITE(3.63)X 

__  63  FORMAT  (  1H  .'I*'J2L 

„(  iw  float  .  . 

0-0.00056 
R-0. 00001 


E9. 1) 


R-R-0.9S 
GO  TO  90 


26 .  80  H1-H1*R*H 

2 

2 

2 

2 

2 

3 

7. 

6. 

9. 

0-0-6.95 

R-R-0.95 

90  IF( (HI  .LE.  0.000001)  .OR.  (HI  .GE.  (0.5-H)))  GO  TO  92 

0. 

1 . 

2. 

N1-N 

Z-INTGRL(H1.R1,R2.N1) 

z-a.o-z 

3 

3 

3 

3. 

N1-1.0 

GAMA  «INTGRL(.H1.R1.R2.N1) 

GAMA-2. 0*A1-GAMA 

3 

3 

3 

CHK»0ABS(E*Bi"?H*H-2 ,0*H*H1 F? .O-GAMaT" 

RHO-((2.0-Z-B2)/CHK)  —  (N/( 1 .0-N)) 
ALPHA-E-B1-(H1--3.0*(H-H1)--3.0)/3.0 

3 
• 

4 

o! 

i . 

N1-0.S 

BETA-INTGHL(H1.R1.R2.N1) 

BETA-2. 0*A1 -BETA 

42. 

43. 

44. 

ni-n7‘(n*i.6)' 

NETA-INTGRL(H1 .R1.R2.N1) 

NETA-2 .0-B2-NETA 

45. 

46. 

47. 

RHIN-1.0/RH0 

M2-RHIN-ALPHA*RHIN*BETA*RHIN--( 1 .0/N)-NETA 

U-0ABS((M1-M2)/M1) 

4 

4 

5 

8. 

9. 

0. 

iP((U  .11.  o.il  .AND.  nr  .GE.  0.6)1  GO  to  60 

WRITE (3, 34)  M2. HI. Ml 

34  FORMAT ( 1H  ,'M2-'.E15.7,'H1«',E15.7.'M1-*.E15.7) 

9 

8 

5' 

5 

5 

S' 

1 . 

2. 

3. 

.... 

5. 

6. 

IP ( Ml  .GT.  M2)  GO  TO  70 

GO  TO  90 

60  KSITA-RH1N-- ( 1 . 0/N)-NETA 

MM- (M2 -KSi TaITT ALPHA  *B  E  TAJ 

V015»( 1.0-DC0S(L-MM))/mm 

HOIS-L*(DSIN( L*MM) )/MM 

— 

5' 

S 

S' 

7. 

9. 

9. 

WRITE (3.291 

29  FORMA T( 1H  . 5X . 'HI '. 10X . 'RHO SX . 'MOMENT '. 9X . 'VOIS 7X . 'HOI S ' ) 
WRITE(3.19)  HI .RHO. M2. VOIS. HOIS 

ii  rOPMATif  1M  " ".'sits  .a  j 

so  to  ao 

ta  WRITEO.0S) . . 

ERROP  MESSAGE 


•TToft»SfiT«rT7i^’’OiTS’T5‘'oSiSUTTiSLE''f6R"DCSlRiEO"STS-ueTUR'fs'1J 
30  CONTINUE 
STOP 


FUNCTION  USED  TO  CALL  INTEGRATION  ROUTINE  DCADRE 


DOUBLE  PRECISION  FUNCTION  INTGRL(H1 . R 1 . R2  .N1  ) 

INTEGER  IER 

"REAL'iDCADRET?".  A.¥’.’AERr7RERR/ERR0R' HT.Cr;C2”C3.'CA.N2. 
1P1.R2.P,R3.RA,H2.N1 

_ COMMON  R3.R4 .N2.P.H2  _  _  _ _ 

EXTERNAL  "f'"" . . . . .  . . 

N2-N1 

R3-R1 


R4-R2 

H2-H1 

RERR-O. 000000001 


A ERR -0.000000001 
A -H2 +0.000000001 
B-R3-R4+H2-0. 000000001 


P-1.5 

CI-OCAORE(F.A.B.AERR.RERR.ERROR.IER) 

A-B+0. 00000001.1 _ 

B-R3+R4*H2 -0.000000001 

C2-DCA0RE (F. A. B.AERP. REPP. ERROR. IER) 


C3-DCA0RE( F . A . B . AERR . RERR . ERROR .IER) 

p-iTs .  . . 

A-B+O.OOOOOOOl 1 

B-2»R3*H2-0.000000001  _ 

C4-DCA0RE  7  F,  A.B .AERR. RERP  .ERROR  .'IER? 

INTGRL-C1+C2-C3+C4 

RETURN 


DOUBLE "PRECI SI ON ' FUNCT I ON  F  <  Y) 

REAL'S  R3.R4.N2.UN0ER1. UN0ER2. H2.Y.P 
COMMON  R3.R4.N2.P.H2 

TfTp  .'it 7 ”*1 .0)  GOfcTTI 

UNDER 1-( 2*R3*Y-Y-Y-H2*H2”2*H2*R3+2-Y-H2)*-0.5 
_ F-(Y— <  1  .0/N2)  )»UN0ER1  _ 

So  fo'i'f .  .  ’ 

14  UNDER2-(R4*R4*2'R3* Y*Y»Y-H2,H2-2-H2*R3+2*H2,Y-R3*R3 )* *0. 5 

F-r  —  (  1  ,0/N2)*UN0EP2  _ _ 

17  REtURN . .  . 

END 

//6.SVSIN  00  •  _ 

0.2  6". 5  6 . 0 156  0': 525 '  0 '.47S  '2 . o  'ao'.O "§075“~07iK02  3  .OF6'7 


PROGRAM vL EARN 


1.  //  JOB  .REGION- 1024K 

2.  //  EXEC  ftvcg 

3  .  .  /7c  .  Si  i IN  DD  •' . . 

4.  c 

8.  C  . . 

. C .  •  SUBROUTINE  FOR  PRINTING  RESULTS  • 

7.  C  - . . . 

a.  .  c 

*  I .  SUBRdUTiNE  RESULT(  Si  rX3 .  V3 ,  RHO,  M) . 

10.  DIMENSION  XORD( 25 ) , YORO( 25 ) 

11.  COMMON  XORO.VORD.XS, V5 

‘  1 2L . ANG-180.6/3 . 1415927-SI . 

13.  X5-V3 

14.  V5--X3 

"  is .  WrIte  i  3  ."asl  "  ang  ,"xs ,vs .rho.  m""v .  . 

16.  *5  FORMAT! 1H  ,*£N0  ROTATION- ‘ , E 12 .4 ,4X. 'XCOR- ' . E 12 . 4 . 4X. 

.17.  _ 1  *  YCOR- * . E 12 . 4 ,4X , 'RHO- ' .E12.4.4X, 'M-  ' ,12) 

18.  . .  MN-M+  i . . . . 

19.  X0RD(MN)-X5 

30.  VORO(MN)-VS 

"  2 1 ". . . RETURN  .  . . .  .  . 

22.  END 

23.  C 

24"  .  C . .  ••••  . 

25.  C  -  SUBROUTINE  FOR  DATA  INTERPOLATION  - 


28.  SUBROUTINE  DTINT(RM1 . ITVPE) 

_ _ REAL  HI. HOIS.  LN 

30.  bi MENS i ON  A( 26.7  V; B ( 2 6 . 7  y.ABTio. 7  j  .  . 

31-  COMMON  /BLK1/RHO.H1.HOIS.VOIS.LN.R1.MM.A.B 

... 32 ...  __  IF(ITYPE.EO.I)  GO  TO  91 

33  . .  iFl  itYPE  mO  )  GO  TO  45  . *. . . . 

34.  45  00  98  1-1. MM 

„  35.  00  98  J-1,7 

36. '“  .  AB  ( 1 ,  J )  -I  ( I J ) . . . 

37.  98  CONTINUE 

_  38 .  __  GO  TO  46 

39.  . .  .  ‘  ‘91  66  99  "I  -1  .  MM .  . . .  .  . .  . 

40.  DO  99  J-1.7 

...41.  AB(  X  .  J)-A<  I ,  J) 

42.  . 99  CONTINUE  . . . . . . 

43.  46  I F ( RM 1  .EO.  0.0)  GO  TO  41 

44.  DO  16  1-1. MM 

45 . . I F(AB(i.lj.GE.RM1)Gb  TO'  tjf .  .  . 

46.  16  CONTINUE 

47.  WRITE (3. 18) 

48.  . 18  'FORMAT!  IH'.' 'THE 'ELEMENT  C  ANNOT  SUBST  A I  NOE  SI  R  ED  MOMENT  ’  j . 

49.  GO  TO  93 

50.  92  K-l-1 

'“'si .  Rif  -  i  rm'i  •  is  t  k  .  iV T/i  abcxTiT-abikTi  ) ) . 

52.  RHC-RIF* ( AB( I . 2 )-AB( K , 2 ) )**E ( K . 2 ) 

53.  H1-RIF-(A8(1 . 3 ) ■ AB( K , 3  )  )*AB ( K . 3  ) 

54.  HblS-RIF*( AB( I . 4 ) -AB i K , 4 ) )*AB(K,4) . 

55.  VOIS-RIF- ( AB< 1 .5 )~AB(K . 5  )  )-AB(K , 5  ) 

56.  LN-ABt 1.6) 

57.  . Ri«AB(  1.7) . . . .  . 

58.  GO  TO  93 

59.  41  RHO-O.O 


Ht-0.0 
H0IS*0.0 
VOIS-O.O 
LN*AB< 1.6) 


93  RETURN 
. ENO . 


•  MAIN  PROGRAM  STARTS  NOW  • 


REAL  LN.RHO.X2. T2.V0IS.H0IS.SI .ANG.R1 .R3.H1 .H3.X3. 
DIMENSION  XORO( 25 ) , YORO( 25 ) ,  A ( 20. 7 ) , B( 20 , 7 ) , IT YPE  1 
1IPLT 1 ( 20) ,N1 ( 20  > . I MOT  1(20) ,RMOM( 20) . RM0M1 < 20) . D< 3 . 
COMMON  X0RD.t0R0.X5.V5 

COMMON  /BLK1/RH0.H1 . HO  I S . VOI S . LN . R 1 . MM . A .B 
REA0d.11)  NN.MM.NU.RL.W1  ,H 

'i  1  FORMAT  ( 31 2 . 3E 1 1 . 4  )  . . 

REA0d.104)  XI.YI.XF.VF 
104  FORMAT ( 4E 1 1 . 4 ) 

read( i.  io)  ( itvpeTuT. iPi(iT.iPLfi(i).Ni ( ij.iMOT i( 

10  FORMAT (51 2) 

REA0(1.12)I (A(1.J).U«1.7),I>1.MM) 

"ia  FORMAT ( 7E1C.3 )  . . . . . 

READ( 1.55)  N2.M2.N7.M7.MX 
55  FORMAT ( 514 ) 

. N2-NN-3 . .  . 

00  26  II-1.NZ 
RMOM( XX )»0.0 


26  CONTINUE 
NZ-NZ-M 


.NSHIF 
) .  IP  1  ( 20) 

•  .4. 13) 


.  1*1 ,NN) 


THE  DATA  FILE  iZ  READ  INTO  THE  DATA  MATRIX 


READ  ALL  THE  ’PHI' 

00  319  1-1.N2 
REA0(40.318)  DU. 1.1. 1,1) 

316  FORMAT ( E 1 1 • 4  ) 

IF(0(I.1.1.1.1).E0  1000.0)  GO  TO  320 
319  CONTINUE 


RE AO  IN  ALL  THE  'RAO' 

320  DO  321  1 • 1 . N2 

1F(0( 1 . 1 . 1 . 1 . 1  )  EO. 1000.0)  GO  TO  324 
oo  322  u«:.m; 

READ)  40 . 323  )  DU.  J.  1.1.1) 

323  FORMAT (£11 . 4  J 

IF(0(l.U.1.1.i».E0  1000.0)  GO  TO  321 
322  CONTINUE 

321  CONTINUE 
C 

C  READ  ALL  the  'THETA 

324  00  325  I-1.N2 

IFIDU.I.l.l.O.EO  1000.0)  CO  TO  330 
00  326  0*2. M2 

1F(0( I. J. 1. 1 . i >.E0. 1000.0)  GO  TO  325 
00  327  K-2.N7 
RE  A0(  40. 328  I  DU. d. *.  1.1) 


120. 

328 

FORMAT (£11.4)  . 

121. 

1MD(  I  .  J.K,  1 .  1 )  .E0  1000.0  )  60  TO  326 

122. 

327 

CONTINUE 

123. 

326 

CONTINUE  . 

124. 

325 

CONTINUE 

125. 

C 

126. 

C 

READ  ALL  THE  'THETA' 

127. 

c 

126. 

330 

DO  331  1-1. N2 

129. 

ZFC0( X . 1.1.1. DEO. 1000.0)  GO  TO  338 

130. 

00  332  U-2.M2 

131 . 

IF(0U,d.  1.  1.1).  EO.  1000  0)  GO  TO  331 

132. 

DO  333  K-2.N7 

133. 

IMD(I.«J.K,1.1).EQ.  1000 . 0)  GO  TO  332 

134. 

00  334  L-2.M7 

135. 

RE AD( 40, 335 )  Od.d.K.L.I) 

136. 

335 

FORMAT! El  1.4) 

137. 

IF(D(I.d.K.L.1).F0. 1000.0)  GO  TO  333 

136. 

334 

CONTINUE 

139. 

333 

CONTINUE 

140. 

332 

CONTINUE 

141  . 

331 

CONTINUE 

142. 

c 

143. 

c 

READ  ALL  THE  'MOM' 

144  . 

c 

145. 

338 

00  339  1-1. N2 

146. 

IF  (0(1. 1.1. 1.1). CO. 1000.0)  GO  TO  348 

147. 

00  34  1  U-2.M2  . . . . 

148. 

IF  (0(1. d.I.I.O.EQ.  1000.0)  GO  TO  339 

149. 

00  342  K-2.N7 

150. 

iF(0(i.d.k.i.i).Eb. iooo.o)  go  to  Hi  . 

151. 

DO  343  L-2.M7 

152. 

IF(D(I.d.K.L.D.EO.  1000  0)  GO  TO  342 

153. 

00  344  M-2.MX  “  . .  . . . 

154. 

RE AD( 40, 345 )  O(X.d.K.L.M) 

155. 

345 

FORMAT ( E 1 1 . 4 } 

156. 

iF(b(i.J.k.L.M).E6. ioob.b)  go  to  343 

157. 

344 

CONTINUE 

158. 

343 

CONTINUE 

159. 

342 

CONTINUE 

160. 

341 

CONTINUE 

161. 

339 

CONTINUE 

162. 

c 

163. 

c 

THE  DATAFILE  HAS  BEEN  READ  IN 

164. 

c 

165. 

348 

RE VINO  40 

166. 

P-0.95 

167. 

E-0.95 

168. 

F-0.95  . 

169. 

G-0.9S 

170. 

c 

171  . 

c 

TO  EVALUATE  MCUR  REOUIREO  TO  GET  THE  LOAD  FROM  INITIAL 

172  . 

c 

POSITION  TO  THE  FINAl  DESIRED  POSITION 

173. 

c 

174  . 

c 

MAIN  LOOP  BEGINS  NO to 

175. 

c 

176. 

c 

177  . 

c 

MAXIMUM  OF  25  ITERATIONS  ARC  ALLOWED 

178. 

c 

179. 

00  56  KK-1,25 

j,YyVjyv%  ■■«>  . >  v  •  AY-  ,s  .  *  /■  .v  /* ■  • 1> ,  •  w>  v  >  .  -  . 


w\ 


v  i»"* >*•  v’  Yr  W  . "v  #"*  i 


j  F (KK . ME . 1)  SO  TO  101 


>s 


r.v 


ft 

§ 


I 


I>> 


ft 


rfij 

VV-J 

ifej 


K3 


Kii 


iso. 

Ml  . 


182. 

C 

THE  HOME NT  IS  APPLIED  ON  THE  LAST  THREE  ELEMENTS  TO  INITIALLY 

183. 

C 

COIL  THE  ARM  AROUND  THE  LOAD 

184  . 

00  102  II-NZ.NN 

185. 

ie»ii*i 

t«6 

RMOM(il )-D( 1.2.2.2.10) 

187  . 

102 

CONTINUE 

188. 

c 

189. 

101 

SI-0.0 

190. 

X2-0  0 

191  . 

Y2-0.0 

M2. 

M-0 

193. 

R3-00 

194. 

H3-0.0 

195. 

Y4»0 . 0 

196. 

IPCK-0 

197. 

1NCK-0 

198 . 

MR-0 

199. 

NSHIF-0.0 

200. 

XORDl  0-0.0 

201. 

YORDl  0*0.0  ' 

202. 

c 

203. 

c 

204  . 

c 

LOOP  FOR  FINDING  THE  END  COORDINATES  AND  ROTATIONS  OF 

205. 

c 

VARIOUS  ELEMENTS 

206. 

c 

207. 

00  20  I-1.NN 

208. 

ITVPE-ITVPEMI) 

209. 

RMI-PMOM(I) 

210. 

CALL  bTiNf'(RMi.iTY>E  j 

211  . 

IF(RMI.EQ.O.O)  00  TO  86 

212. 

IP-IP1 ( I ) 

213. 

CO  TO  87 

214. 

86 

IP--2 

215. 

87 

IPLT-1PLTKI) 

216. 

N-NKI  ) 

217. 

I MOT - 1 MOT  1(1  ) 

218. 

IF<! .EO. 1)  GO  TO  17 

219. 

IFUDIR.EO  IMOT)  CO  TO  17 

220. 

IF1SI .LT.O.O)  GC  TO  19 

221  . 

HR -MR-  1 

222. 

c 

THE  NUETRAL  AXIS  IS  SHIFTEO  ACCORDINGLY 

223. 

c 

224  . 

IF(ICIR.EO.O)  CD  TO  21 

225. 

IF(MR.EO.I)  NSHIF-R3-H3*R1*H1 

226. 

X2*X2*NSHIF*S1N(  SI  I 

227. 

Y2*Y2*NSHIF*C0S(SI  ) 

228 

NSHIF-0.0 

229. 

c 

230. 

c 

TO  CHECK  IF  PRESSURIZED  OR  NOT 

231  .■ 

.  c 

232  . 

IMIP.EO  -2)  GO  TO  23 

233. 

X2-X2-1 LN-HDIS ) *COS ! SI )-VDIS*SIN( SI ) 

234. 

V2- Y2- ( LN-HOI S ) *S1N( S I )*VDI S*COS(  SI ) 

235. 

SI-SI-LN/RMO 

236 

IF (  S3 .LT.O.O)  GO  TO  24 

237. 

X3-X2-(R1*H1 )*SIN(SI  ) 

238. 

Y3-Y2-1R1-H1  )*CCS< SI ) 

239. 

GO  TO  25 

Hi 


i 


a 


I 


X3»X2*(»i*H1 )»SIN(-SI ) 
V3»Y2- ( R 1*H1 )*CDS( -51  ) 

M«M*1 

CALL  RESULT (SI. X3 . V3 . RHO .Ml 

R3-R1 

H3*H1 

60  TO  20  . 

X2-X2*LN-C0S(SI  ) 
Y2«Y2-LN»$IN(S1) 

X3«X2-(R3*H3)»SlN(Si ) . 

Y3»Y2-(R3*H3)»C0S(SI ) 

M«M*1 

CALL  RESULT (SI .X3.Y3 .RHO. M ) 
60  TO  20 


WHEN  REVERSE  CURVATURE  MOTION  IS  IN  OPP  DIRECTION 


IF(MR.EO.I)  NSHI F»R3*H3*R 1*H1 

X2 ■ X2+NSHI F  •  S I  Nt  S I)  . 

V2»V2-NSHIF-C0S(S1 ) 

NSHI F -00  _ 

IF( IP. 16.-2)  60  TO  33 . 

X2»X2+(LN-HDIS)*C0S(SI )*VDIS*SIN(SI ) 

V2*V2*( LN-HD1 S )*5IN( SI )-VDIS*COS( SI ) 

SI »  SI - LN/ RHO 

IF(SI  LT .0.0)  60  TO  34 

X3»X2* ( R 1*H1 )*SINI SI ) 

'"V3«V2*( R i *H i  )'*cos<  si) 

60  TO  3S 

X3«X2*(R1*H1 )*SIN(  -SI ) 

V3« Y2*( R 1*H1 ) *COS(  -Si ) . . 

M»M*1 

CALL  RESULT ( SI . X3 . Y3 , RHO ,M ) 

R3-R1 . . 

H3-H1 

60  TO  20 

X2»X2-LN*C0S(Si) 

Y2* Y2*LN*S1N( SI ) 

X3«X2- ( R3-H3 ) *SINf  SI ) 

Y3*V2*(R3*H3 )*COS(  SI ) . 

W»M—  1 

CALL  RESULT(SI .X3.V3.RH0.M) 

60  TO  20 
S1*-SI 

IOIR-IMOT 
00  t 5  U-i.N' 


■  'HI  11 


J uujjmiuiu j u. imii imj.jp >u mj  j i .uipjmi uuuuw.< v \ 9 juuim 


300. 

X2»X2-NSHIF-SIN( S3 ) 

301  . 

Y2“ Y2*NSHI F *C0S( SI ) 

302 

39 

X2«X2*< UN-HOIS )»C0S< SI )-VDIS*SIN(SI) 

303 

Y2»V2- ( LN-H01 S )*S1N( 51 ) - VD1S*C0S( SI ) 

304. 

SI«51*LN/RH0 

305 . 

X3»X2*<R1*H1 )»SIN(SI ) 

306. 

Y3»V2*(R 1*H1 )*C0S<S1 )->4 

307. 

R3»R  1 

308. 

H3-H1 

309. 

CALL  RESUL7(Si.X3.Y3.RH0.M) 

310. 

GO  TO  IS 

311 . 

40 

M-M*1 

312. 

i  r '( iNC'it .  Ed .  6)  Y4— ( R  1*H1 ) 

313. 

INCKb 1 

314. 

IF  (M  .EO.  1)  GO  TO  49 

315. 

NSMiF»(R1*M1)-(R3*H3) 

316. 

X2bX2-NSHIF«S!N(SI ) 

317. 

Y2«V2*NSHIF*C0S(SI ) 

318. 

49 

X2*X2*( LN-HDI S ) *COS( SI )>V01S*SIN(Si ) 

319. 

Y2-Y2*(LN-H0IS)*S1N<SI I+VOI S*C0S( SI ) 

320. 

SJ-SI*LN/RHO 

321  . 

X3-X2*(R1*M1 y«siN(si ) 

322. 

Y3«Y2-(R1*H1)»C0S(SI )-Y4 

323. 

R3“R1 

324. 

H3-H1 

325. 

CALL  RESULT ( SI , X3 . Y3 , RHO ,M) 

326. 

GO  TO  15 

327. 

C 

328. 

C 

TO  CHECK  OIRECTION  OF  MOTION  OF  STRUCTURE 

329. 

c 

330. 

30 

i'F( iPCk  .EO.  0)  GO  TO  S3 

331. 

JF(IMOT  .EO.  -3)  GO  TO  51 

332. 

c 

333. 

c 

THE  TRANS FORK AT 1  ON  FOR  DOWNWARD  TRANSLATION 

334. 

c 

335- 

X2«X2-»LN*C0S{SI  ) 

336 

Y2»Y2-LN*SIN(SI) 

337. 

X3«X2*( R3*H3 )*SIN(SI ) 

338. 

Y3*V2*( R3*H3 ) *COS( SI )• Y4 

339. 

M>M*1 

340. 

INCK* 1 

341  . 

CALL  RESULT ( SI . X3 . Y3 . RHO , W ) 

342. 

60  TO  15 

343. 

51 

X2bX2*LN*C0S(S1 ) 

344  . 

Y2bY2*LN»SIN(SI  ) 

345. 

X3«X2*(R3*H3 )*S!N(  SI  ) 

346. 

Y3bY2-(R3*H3)*C0S( SI )*V4 

347. 

MbN*  1 

348. 

lNCK* 1 

349. 

CALL  RESULT(SI.X3.Y3.RhO.K) 

350. 

GO  TO  15 

351  . 

c 

352. 

c 

WHEN  elements  IN  THE  BEGINNING  are  not  pressur 

353. 

c 

354. 

53 

X2«X2*LN 

355. 

Y2-0.0 

356. 

MbM*1 

357. 

INCKb  1 

358. 

X3-X2 

359. 

Y3»Y  2 
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CALL  RE  SUIT ( SI , X3 , Y3 . RHO .  M) 

IS  CONTINUE 

LOOP  POR  IDENTICAL  ELEMENTS  IS  COMPLETE 
20  CONTINUE 

LOOP  FOR  FINDING  END  COORDINATES  FOR  ALL  THE  ELEMENTS 
IS  COMPLETE 


MN»M* 1 

CALL  PICS! Z( 20.0. 1 1 .0) 

CALL  0RIG1N1 5 .6.  10.5)  . . . . . . 

CALL  FACT0R(2. 0.2.0) 

CALL  AXIS(-3. 0.0. 0.230. 'POSITION  OF  STRUCTURE  IN  XDIRECTION 

'1.6. 6. 2.0. -3. 6. 20.0) 

CALL  AX1S( -3. 0.0. 0.20.0. 'POSITION  OF  STRUCTURE  IN  YOIRECTION 
1.-90. 0.2. 0.0.0. -20.0) 

"CALL  symbol (2. 6. -4. 0.6. 2.  . . 

1 'DEFORMED  ROBOTIC  ARM  ON  INTERNAL  PRESSURIZATION  '.00) 

CALL  PLOT(XORD.YORD.MN) 


WRITE  (  3  ,  B4  )  (  RMOM (  X'Y  )  .  I Y *  1  , NN  ) 

B4  f 0RMAT11H  .* MOMENT  AT  PRESENT* • .E12.5) 


NOW  TO  CHECK  THE  CURRENT  POSITION  OF  THE  LOAD  W.R  T  THE 
THE  FINAL  DESIRED  POSITION  AND  THEREBY  VARY  THE  MOMENTS 
IN  THE  ELEMENTS 


IF(KK.EO.I)  GO  TO  69 
P1»P 
"E1-E 
F1-F 
G1-G 


TO  EVALUATE  THE  COORDINATES  OF  THE  CENTER  OF  THE  LOAD 


. XB-XS-  ("RLiR  1*H j -COS  ( SI )  '• . . 

VB-Y5*(RL*R1*H)*SIN(SI ) 

WRITEO.85)  XB.yB.xF.vf 

85  FORMAT (  1H  .  'XE*'  .E12.5.5X.''YB»'''.Ei'2'.5.5X"'';XF»;'.E10.3.5X. 

1 '  VF* • . E 10 . 3 ) 

XD1S-XB-XF 

YDiS-YB-YF 

IF ( 1  ABSt XDIS).LE.1.0). AND .(ABS(YDIS).LE.I.O))  GO  TO  60 

NOW  TO  CHECK  THE  CURRENT  POST  ION " OF  LOAD  W.R.T  FINAL  POSITION 

IF((XDIS.LT. 0.0). AND.IYOIS.LT. 0.01)  GO  TO  62 
IFUXCIS.LT. 0  0). AND  1YDIS.GT  .0.0)  )  GO  TO  63 
IF1IXD1S  GT.O. OJ.AND.lYDIS.LT. 0.0))  GO  TO  64 


WHEN  LOAD  IS  TO  THE  RIGHT  A NO  HIGHER 


NA»NU* 1 
NE-NN-3 

P1»Pl/12.0*1 1NN-3I-NU)) 
P2»P1 

00  65  i-NA.NB 

RM0M1 I )«RM0M11 )*RM0M( I )»01 

Pl»Pl*P2 


430 

421 

422 

423 

424 

425 

426 

427 

428 
"429 

430 

431 
"432 

433 

434 


435.  GO  TO  83 

436.  C 

437.  C  WHEN  THE  LOAD  IS  TO  THE  LEFT  AND  HIGHER 

. 438.  ™C . 

439.  €3  F 1*F 1/(2 .0*NU) 

440.  F2-F1 

. 441 ; .  DO  67  T-T.  NU . .  * . . 

442.  RMOM(I)>RMOM(1)-RMOM(1)*F1 

443.  F 1«F 1+F2 

444" . 67  CONTINUE . . . . . 

445.  F-F-F»0.05 

446.  _  GO  TO  83 

"  ’447  .' . C . . . . . . .  . 

448.  C  WHEN  LOAD  IS  TO  THE  RIGHT  AND  LOWER 

_449.  c 

'450. . . . 64"  CT-Gl/(  2"'6-NU) . ~ . 

451 .  G2»G1 

452.  00  68  I-1.NU 

453.  ' .  RMOM(  I  )»RMOH(  I  )*RwbM(  i  )*(ii . . . 

454.  G1«G1*G2 

455.  68  CONTINUE 

'  456. . G»G-G*0 . 05 . . . . 

457.  GO  TO  83 

458.  C 

459.  . C . THE  DATA  MATRIX  is  SEARCHED  TO  PROVIDE  THE  GUESS  VALUE . 

460.  C 

461.  C  TO  EVALUATE  PHI .RAD . THETA  6  DIS  FOR  THE  GIVEN  COORDINATES 

462.  . 69  RAD»(ABS«XI )*'2.0«ABS{vi >»*2  b)-*b.5 

463.  D1S*( (AES(XF-XI  )  )  »  *2 . 0*1AB5(VF-YI ) )**2.0)**0.5 

464.  THETA»ACOS< (XF-Xl )/DIS) 

465.  IF(xi .LE.O.O)  GO  TO  360 

466.  PH1«ATAN( X1/ABS1 YI  )  ) 

467.  GO  TO  361 

468.  .  360  PHI ■-ATAN(Xi/Vl ) 

469.  C  THE  SEARCHING  STARTS  NOW 

470.  C 

471.  C  TO  FIND  THE  CLOSEST  PHI 

472.  361  1«0 

473.  362  I-I-1 

474.  IFCPH1.GT.DU. 1.1. 1.0)  GO  TO  362 

475.  K-l-1 

476.  DIFFaABSC  PHI -DCK.1.1,1.1)) 

477.  D!FF2*ABS(D(I . 1. 1. 1. 1 )-PHI ) 

478.  IFC01FF .G7 .DIFF2)  GO  TO  363 

479.  I-K 


65  CONTINUE 
P-P-P»0.05 
GO  TO  83 

C 

C  WHEN  LOAO  IS  TO  THE  LEFT  AND  LOWE* 
C 

62  NC-NU* 1 

El»E1/((2.0«(NN-3)-NU)) 

E2-E1 

DO  66  1-NC.NO 
RMOMC I )>RMOM( 1 )-RMOM( I ) - E 1 
.  E  1“£  i*E2 .  . . . 

66  CONTINUE 
E*E-E*0 . 05 


'  4S0 . 

481  . 

482  . 


4es. 

486  . 

487. 

488. 
4897 
490. 

491  . 

492  . 
492. 
494  . 

495. 

496. 

497. 
498'. 

499. 

500. 

501  . 

502  . 
503. 
504 

505. 

506. 

507. 

508. 

509. 

510. 
511  . 

512. 

513. ' 

514. 
$15. 

516. 

517. 

518. 

515. 
520. 
521  . 

522. 

523. 

524. 

525. 

526. 

527. 
528  . 

529. 

530. 
531  . 

532. 

533. 
534  . 

535. 

536. 

537  . 

538 
539. 


C 

C  TO  FIND  THf  CLOSEST  BAD 

363  U«1 

. 364  J-J*1 . 

1F(RAD.GT.D(I.J. 1,1.1))  GO  TO  364 
L  *U*  1 

I F ( (RAD-D( I .L. 1. 1 . 1 ) j  GT .(D(I. j. 1.1. 1 )-RAD j )  GO  TO  385  . 

U-L 

C 

C  to  FIND  THE  CLOSEST  THETA 

365  K» 1 

366  K*K*1 

. IF ( THETA .6T.D< I .d.K.  1. 1 ) )  GO  TO  366 

M»K-  1 

1 F ( ( THETA”D( I.J.M.1,1)), GT ,(D(1. d.K. 1,1 )-TH£T A ) )  GO  TO  367 

K-M  ' 

C 

C  TO  FIND  THE  CLOSEST  DIS 

. 367  N-1 . 

366  N»N*1 

I F  ( D I S  .  GT  .  D ( I  .  J ,  K  ,  N ,  1  )  )  GO  TO  368 
L-N-1 

IF((DIS-DO.d.K.L.1)).GT  (DO.J.K.N.1)-DIS))  GO  TO  369 
N-L 

C  . 

C  TO  TRANSFER  THE  GUESS  VALUE  ' MOMENTS '  IN  THE  MEMBERS  TO 
C  START  THE  ITERATIVE  PROCEDURE 

"c . 

369  DO  370  M1-1.NN 
M«M1*1 

. . ~RM3M(  Ml  )  -D(Xd".  K  ,  N .  M  ) . . . . 

370  CONTINUE 
C 

C . THE  "GUESS" MOMENT  VALUES  HAVE  BEEN  LOCATED . 

C 

83  VRITE(3.57)kk 

.  57  FORMAT ( 1H  .'ITERATION  NO1 .12.'  IS  COMPLETE ') . 

56  CONTINUE 

C  MAIN  LOOP  ENDS  NOV 

"C  the  MCUR  HA 5  BEEN  EVALUAtED  VI thin  25  ITERATIONS 
C 

c 

. WRITE (3. 97) 

67  F0RMAT(1H  ,'NOT  POSSIBLE  TO  REACH  GOAL  IN  25  STEPS') 

GO  TO  172 

C '  " " 

60  CALL  PICSIZ(0.0,0  0) 

C 

C*' . . . . . . . . 

C 

C  . . . . . . . . . . 

C  •  THE  DAtABASE  IS  CHECKED  AND  THE  PRESENT  RESULTS  ARE 

C  •  INSERTED  IN  THE  PROPER  PLACE  IN  THE  DATABASE  IF  NEEDED  • 

C  — . . . . . . . . . 

C 

DIS-((xe-XI)**20*<vE-vi)-*2.0>**0  5 
theta»acos((xb-xi )/d:s> 

C  TO  PUT  PHI  IN  PROPER  PLACE  IF  REOUIRED 

I»0 

127  1-1*1 
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TV*  . 


rv  .j"  rv 


540.  IF(A8S<PHI-D(I.1.1.1.1)).LE. 0.00009)  SO  TO  126 

541.  C  NO  NEED  TO  INSERT 

542  1F(PHI.GT.D(I.1.1.1.1))  GO  TO  127 

543.  "  "C  "TO  INSERT  AND  PUSH  THE  DAf A  DOWN  SUCCESSIVELY 

544.  IP-I 

545.  00  150  I-1.N2 

546.  . 1F(D(I . 1. 1. 1. 1)  £0. 1000  6)  GO  TO  151 

547.  150  CONTINUE 

548.  151  1K>I*1 

"549.  . . .  IE (  ik.GT  1N2  )  60  TO  169 

550.  C  THE  INSERTION  AND  SHIFTING  LOOP  STARTS  NOW 

551.  0(IK.1.1.1,1)-D(I.1.1.1.1) 

552.  IK»ik-1 

553.  1-2-1 

554.  IP5-IP-1 

MS. . 167  iF(i.C0.iP5)  GO  TO  166 . 

556.  0(IK.1.1.1.1)»0(I.1.1.1.1) 

557.  DO  152  J«2.H2 

558.  . oiXK.J.  i.i.  i)-6(i.j.  i.T.  ij . .  . . 

559.  IF(D(I.J.1.1.1).E0. 1000.0)  GO  TO  153 

560.  152  CONTINUE 

"  56 1". .  1 53  dd-d-  1 . . . . . . . . . 

562.  DO  154  U-2.UJ 

563.  00  155  K*2.N7 

'564'.  . 0(  ik.d.K.  1 .  1  ) - 0 < i  , U.K,  1 ,  i  r~ . 

565.  lF(0(l.d.K.1.1).E0. 1000.0)  GO  TO  156 

566.  155  CONTINUE 

567.  . 154  CONTINUE  . . .  . . . . 

568.  156  KK-K-1 

_  569 „  00  157  d«2.dd 

570.  '  oo“iSI"k-a.kK . .  .  . 

571.  00  159  L-2.M7 

. 572. __  0UK.J.K.L.1)*0a.0.K,l.1) 

573 .  1 F  ( D<  I .  d .  k .  L .  i ) .  60  1000 . 6)  GO  TO'  160 ~ . . . . 

574.  159  CONTINUE 

575.  158  CONTINUE 

576.  . 157  CONTINUE  . . 

577.  160  LL-L-1 

...578. .  00  161  d«2.dd 

579.  '  'DO  162  K«2.KK  .  . 

560.  00  163  L-2.LL 

581.  DO  164  LM«2.MX 

"  582.  .  0< IK. d.K.L.LM) >0(1 .d.K.L.LM) 

583.  IF(D(I .d.K.L.LM) .EO  1000.0)  GO  TO  165 

584.  164  CONTINUE 

"  "585. . 163  CONTINUE 

586.  162  CONTINUE 

587.  161  CONTINUE 

588.  . C  TO  GO  UP  ONE  INDEX  AND  PUSH  DOWN 

589.  165  I F ( IK . EO .IP)  GO  TO  166 

590.  lw.lK-1 

591  .  . I  - 1  -  1 

592.  GC  TO  167 

593.  C  SHIFTING  IS  COMPLETE 

594.  166  S( IK. 1 . 1 . 1 . 1 )>PH1 

595.  0( 1K.2. 1 , 1 . 1 )»RAO 

596.  C(IK.3.1.1.l)«lOOC.O 

597.  .  C(lK,2.2,i.l ) "THETa 

598.  S( IK . 2 . 3 , 1 . 1 )• 100C  0 

599.  0( IK. 2.2.2.  O-DIS 
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D(IK , 2 .2 . 3 . 1 )• 1000 . 6 

MC-NN* 1 

00  168  I6«2.MC 

kk-IB-1 

OdK.2.2.2.  IB)-RMOM(KK) 

168  CONTINUE 

. . . . 

MC"MC*1 

D(IK,2,2.2.MC)* 1000.0 

INSERTION  COMPLETE  GO  TO  WRITE  FILE  AREA 

GO  TO  170 

ERROR  MESSAGE  IF  FIRST  DIMENSION  WAS  EXCEEDED 

169  WRITE ( 3, 171 ) 

171  F0RMAT(1H  .'THE  FIRST  DIMENSION  WAS  EXCEEDED') 

GO  TO  END  OF  PROGRAM 

GO  TO  172 

IF  PHI  WAS  NOT  INSERTED  SEARCH  TO  SEE  IF  RAD  NEEDS 

TO  BE  INSERTED  IN  THE  DATABASE 

126  J"1 

174  U"U*1 

iF(ABS(RAb-b(i.J.1.1.i)).LE.b  009)  GO  TO  173 

IF(RAO  GT .0(1. U. 1.1.1))  GO  TO  174 

TO  INSERT  RAO 

JP-U 

00  175  U*2.M2 

IFlDd.U.  1.1.1).  EO.  1000.0) 

GO  TO  176 

175  CONTINUE 

176  UK-U-M 

1 F ( UK . GT . M2 )  GO  TO  177 

d( i.Jk, i . i .i y-o(i.u.i.i.i) 

JK" JK- 1 

JP1-JP-1 

188  IF(U.EO.UPI)  GO  TO  18? 

D(I.UK.1,1.l)"0<I.J. 1,1.1) 

DO  178  K-2.N7 

Dd.JK.K,  1,  1)*0(I.i).K,  1.1) 

IF(Dd.J.K. 1.1). EO.  100C.0)  GO  TO  179 
179  CONTINUE 

179  KK-K-1 

DO  180  K"2.KK 

. 00  181  L«2.M7  . 

Od  .  JK.K.L.  1)*0d.U.K.l.  1 ) 
!F(Dd.J.K,L.1)E0.  1000  0)  GO  TO  182 

181  CONTINUE 

180  CONTINUE 

182  LL-L-1 

00  183  K-2.KK 
00  184  L • 2 . Lt 
00  18S  LM.2.M1 

0(1 .JK.K.L .  LM)"Dd . J.K.L.LM) 

IF<0< I . J.K.L.LM) . EO. 1000  0)  GO  TO  186 
185  CONTINUE 
184  CONTINUE 

183  CONTINUE 

GO  UP  4 NO  PUSH  DOWN  AGAIN 


186  IMUK.EO.JP)  GO  TO  187 
UK  »  JK  -  1 
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J.J-1 

GO  TO  188 

1S7  D(I.JK.1.1.1)-RAD 
D(I.dK.2.i.1)*THETA 
0(1 . JK.3, 1 . 1 >*1000.0 
0d.JK.2.2.1)*DIS 
D(I.UK.2.3.1)*1000  0 
MC*NN* 1 

00  189  18*2. MC 
KK-I8-1 

Dd  .  JK.2.2.IB)-RM0M(KK) 

189  CONTINUE 

MC-MC+1 . . . : . 

0(1  .  JK. 2. 2. **0*1000  0 
GO  TO  170 

C .  the' SECOND  DIMENSION  WAS  EXCEEDED . 

177  WRITEd.190) 

190  FORMAT! 1H  .'THE  SECOND  DIMENSION  HAS  BEEN  EXCEEDED') 

C .  GO  TO  THE  END  OF  PROGRAM  WITHOUT  REWRITING  THE  DATAFILE 

GO  TO  172 


IF  PHI  &  RAD  WERE  NOT  INSERTED.  CHECK  FOR  THETA 
K*1 


K«K-*1 

XF(ABS(TMETA-0(I.J.K.1.1)).LE. 0.00009)  GO  TO  191 
IF (THETA .GT .O(I.U.K.I.I))  GO  TO  192 


TO  INSERT  THETA 
KP-K 

DO  193  K-2.N7 

TFfDdTu.'k. i.i).Ib.‘ib6c  b)  go  to  iST" 

CONTINUE 

KK-K«1 

'  IF(*CK  "<ST  .N7)  GO  TO  195 . . 

D(I.d.KK.1.1)*D(I.U.K.1.l) 

KK-KK-1 

k«k-  i . . . . 

KP1-KP-1 

IF(K.EO.KPl)  GO  TO  201 
Dd.d.KK.I.U-Dd.d.K.I.i) 

DO  196  L-2.M7 

D(I . J.KK.L. 1 )*0( I .U.K.L .  1  ) 

IF(D(I.J.K.L.1).EC. 1000.0)  GO  TO  197 

CONTINUE 

LL-L-1 

00  198  L-2.LL 

DO  199  LM-2.MT 

D( I.U.KK.L.LM1*D( 1 . J.K.L.LM) 

iF(0(I.U.k,L.LM).E0. 1000  0)  GO  TO  200 

CONTINUE 

CONTINUE 

GD  UP  AND  PUSH  DOWN  AGAIN 

IF(KK.EO.KP)  GO  TO  2C i 

K*K-1 

KK*KK-1 

GO  TO  202 

D(  I , U . KK . 1 . 1 ) • T HE T  A 
0( I . J.KK.2, 1 )*DIS 
Dd.J.KK.  3.1)*  1000.0 
MC*NN*1 


we 


721.  M2  * 1 B - 1 

722.  0(1 .J.KK.2.IB)«RMQM(M2) 

‘723.  .  272  CONTINUE 

724.  MC«MC«1 

725.  D( I . J.KK . 2 . MC )* 1000 . 0 

726.  C  . INSERTION  COMPLETE  00  TO  WRITE  FILE  ARE* 

727.  GO  TO  170 

728.  C  THE  THIRD  DIMENSION  HAS  BEEN  EXCEEDED 

'729'.  . 19*  WRITE (3. 203) 

730.  203  FORMATOH  .'THE  THIRD  DIMENSION  HAS  BEEN  EXCEEDED*) 

731.  C  60  TO  END  OF  PROGRAM  WITHOUT  REWRITING  THE  DATAFILE 

'  732.  .  GO  fO  172 . 

733.  C 

734.  C  IF  PHI. RAD  6  THETA  WERE  NOT  INSERTED,  CHECK  FOR  DIS 

735  ! . C . 

736.  191  L«1 

737.  205  L*L* 1 

"73*  . iF(  ABS(OiS-D(I .J.K.L.  1  jV.LE. 6.00009)  GO  TO  237  . 

739.  XF(DIS-GT . D( I .J.K.L, 1))  GO  TO  205 

740.  C  TO  INSERT  DIS 

“741. . LP*L  .  .  ”  . "7  . .  . .' . 

742.  00  206  L-2.M7 

743.  IF(0(I.J.K.L.1)-EC. 1000.0)  GO  TO  207 

744. ' . 206  CONTINUE 

745.  207  LK«L*1 

746.  IFCLK.GT.M7)  GO  TC  20B 

"747; . . 0(1  .U.K.LK.  1  )«0(i  .J.K.L.  iT 

74*.  LK-LK-1 

749.  L«L-1  _  ' 

"750" . . LP1-LP-1  . . 

751.  212  IF(L.EO.LPI)  GO  TO  211 

752.  OU.J.K.LK.D-Od.J.K.L.D  __ 

'  753; . ""00  209  LM-2.MX  . . 

754.  0(1. J.K.LK.LM)«0(I .O.K.L.LM) 

t«5.  IF(D(I.J.K.L.LMi.EO. 1000.0)  GO  TO  210 

“'756. . 209  CONTINUE 

757.  C  GO  UP  AND  PUSH  DOWN  AGAIN 

758.  210  IFILK.EO.LP)  GO  TO  211 

759.  .  LK-LK-1 

760.  L-L-1 

761.  60  TO  212 

762  . . .  211  MC*NN*1 

763.  DO  213  IB-2.MC 

764.  M2-IB-1 

765.  0( I . J.K.LK.JB)-RMDM(MZ) 

766.  213  CONTINUE 

767.  MC  •MC'*  1 

768.  "  .  D(  I.  J.k.LK.MC)*  1000  0  . . 

769.  C  INSERTION  COMPLETE  GC  TO  WRITE  FILE  ARF A 

770.  60  TO  170 

7*>i.  C  THE  FOURTH  DIMENSION  WAS  EXCEEDED 

772.  20B  WR1TE(3.214) 

773.  214  PORMATliH  .'THE  FOURTH  DIMENSION  HAS  BEEN  EXCEEDED*) 

774.  C  GO  TO  END  OF  PROGRAM  WITHOUT  REWRITING  DATAFILE 

775.  60  TO  172 

776 .  C 

777.  . C 

77*.  C  TO  OVERWRITE  THE  DATAFILE  WITH  THE  NEW  OAT* 

779.  C 


c 

c 


I 

*4 

$ 

it 

$ 

m 

i 

i 


WRITE  ALL  THE  'PHI'  FIRST 

DO  215  1-1. N2 
WRITE(40.2i€)  Dd.i.i.i.i) 

FORMAT (Ell .4) 

IF(0(I.1.1.1.1)EO  1000  0)  GO  TO  217 
CONTINUE 

WRITE  ALL  THE  'RAO' 


00  2 18  !>1.N2 

IF(0(I.1.1.1.1)E0. 1000.0)  GO  TO  219 

00  220  U*2,M2 . 

WRITE (40, 22 1 )  0(1. J. 1.1.D 
F0RMAT(E11 .4) 

iMOti  .U.  1. 1.1). EO.  lOOO.O)  GO  TO  218 
CONTINUE 

CONTINUE .  . 

WRITE  ALL  THE  'THETA- 


DO  222  I-1.N2 

1FI0(1.1,1.1.1).E0. 100C . 0)  GO  TO  223 
00  224  J«2,M2 

IF(D(I.J.1.1.1).E0.1006.0j  GO  TO  222 

00  225  K-2.N7 

WRITE  (40 . 226  )  Dd.J.K.1.1) 

FORMAT ( Ell . 4 ) 

IFtDd.d.K.  1,  1).E0. 1000.0)  GO  TO  224 
CONTINUE 
‘  CONTINUE 
CONTINUE 


WRITE  ALL  THE  'DIS' 

00  228  I ■ 1 . N2 

IF(0<I.i.1.1.i).EO. 1000.0)  GO  TO  22* 
00  230  U-2.M2 

IF(D(I.d.1.1.1).E0. 10000)  60  TO  228 
DO  231  K-2.N7 

IF(0(I.J,K.1.1).EO. 1000.0)  GO  TO  230 
DO  233  L-2.M7 

WRITE ( 40,234 )  O(i.d.K.L.i) 

FORMAT  (EH.4) 

1F(D(I.0.K,L.1).E0. 1000  0 )  GO  TO  231 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

WRITE  ALL  THE  ELEMENT  MOMENTS 
00  236  I » 1 . N2 

Ir(0(I.1,1.1.1).E0. 1000 .0)  GO  TO  237 
DO  238  J>2.M2 

IF(0(I.J. 1.1,1). CO. 1000.0)  GO  TO  236 
00  239  K*2.N7 

!F(D( I . J.K. 1 . 1  )  EO. 1000  0)  GO  TO  238 
00  241  L-2.M7 

IF(0( I .d.K.L. 1 )  EO  1000. 0)  GO  TO  239 


%  wmsmmmmmmmm. 


•40.  .  00  243  M-2.MX . 

•41.  WRITE ( 40. 244 )  D(I.J.K.L.M) 

•42.  244  F0RMAT(E11 .4) 

•43.  . XF(D(1 .O.K.L.M)  EQ.  1000. 6)  GO  TO  241 . 

844.  243  CONTINUE 

■45.  241  CONTINUE 

846.  239  CONTINUE 

•47.  238  CONTINUE 

•48.  236  CONTINUE 

849 . . c . ; . . 

•SO.  C . . . . . . . 

•SI .  c 

•52  . . C . *••* •»•*•••  -•  * -* •••  *•  •  ••••*••  • 

•S3.  C  *  TO  CALCULATE  THE  ADDITIONAL  MOMENT  DUE  TO  THE  LOAD  * 

•54 .  C  •  TO  CALCULATE  The  X-COORDINATE  OF  POINT  C 

'■"•«.  C . . . ........................................ 

•56.  C 

•57.  237  GAMA*A8S( SI -3*LN/RH0 ) 

'"•58.  RM3Wwi«(RL*R1*H) 

•59.  XC»XB*(RL*R1*H)*C0S(GAMA) 

•60.  RM2*W1‘fRL«Rl+H)»C0S(GAMA) 

. 861 . . 00  70  J-1.NU 

•62 .  RMOM( J).RMOM( J)-RM2*W1*(XC-X0RD( J) ) 

•63.  70  CONTINUE 

"■'•64."  NF»NU*1  . - . 

•6S .  NG.NN-3 

•66.  00  71  J-NF.NC 

. *67  . . RM0M(  J  ).RM0M(  j  )*RM2-W 1  •  (  XC-XOROi  J  )  j . . . 

•68.  71  CONTINUE 

•69.  NH-NN-2 

■■""•70 1 . . .  DO  72  «i*NH.NN  .  .  . . . . 

•71.  RMOM( J)«RMOM( J)*RM3 

•72.  __  72  CONTINUE 

"■""•73. . <a.s»y"i KMSMl'f If7f>  W)  ”  "  .  . . 

•74.  SB  FORMAT! 1H  .‘TOTAL  MOMENT.  '.£10. 3) 

•75.  172  STOP 

"  976 . .  END 


•77. 

/• 

•78. 

//G.FT4CF001 

DO  DSN*DU  006 . 2U2 142  UMESH . LOATA ,D1 SP* 

OLD 

•79. 

//G.DUKEPLOT 

00  STSOUT-C 

■80. 

//G  SYSIN  00 

• 

•81  . 

1010  5  0. 1944E*01 

0.26*01  0.  116-01 

•82. 

-0  *95 i£*0i 

-0. 14E*C2  0.81*01  -0 . 9E*01 

•83. 

1*2-1  1-3 

•84 

1*2-1  1-3 

US 

1*2-1  1-3 

•86. 

1*2-1  1-3 

•87  . 

1*2-1  1-3 

•88. 

i*?»i  1*3 

■89. 

i*2*i  1*3 

■90. 

1*2*1  1*3 

891  . 

1*2*1  7*3 

892. 

1*2*1  1*3 

•93. 

C. 2S*02 

0 

115E*02 

0  5E-C2  0. 101E-01 

0. 173£«0C 

0.2E*Ci 

0 . 5E*00 

994 

0  256*02 

c 

764E*Ci 

C  499E-02  C  2281-01 

0  266*00 

0 . 25*01 

0 . 5E*00 

•95. 

0  36*02 

0 

573E*01 

0  49BE-C2  0  4046-01 

0  346E*00 

0 . 2E»01 

0.5E*00 

■96 

0  356*02 

0 

4581*01 

0  4976-02  0.630E-01 

0  431*00 

0  ■  2E*01 

0. 56*00 

•97. 

C  4|*02 

0 

382E*0i 

0  4966-02  0  90 IE -01 

0.5126*00 

02E»0i 

0. 5E*00 

899 

C  451*02 

0 

327E*0‘ 

0  4956-02  0. 123E*00 

0 . 593E  *00 

0.2E*0i 

0.56*00 

■99. 

0  5E«02 

0 

2ITE*0i 

0  494E-02  0. 158E*00 

0 . C69E*00 

0.26*01 

0  56*00 
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6 . S5E+02  0.255E*01  04931-02  0.199E*00  0.745E*00  0.2f*01  6  5E*00 

0.«E*02  0.229E*01  0.492E-02  0.24SE*00  0.ei9E*00  0  2E*01  0  5E*C-0 

O- 7E*02  0. 191E*01  0  491E-C2  0.346E-0G  0.955E*00  0.2E-01  0.5E*0C 


